Influence of DE-cluster refinement on numerical analysis of rockfall experiments

A numerical analysis is validated against a Swiss Federal Commission for Technology and Innovation (CTI)—frame impact experiment conducted by the Swiss Company Geobrugg. The discrete element method is used to simulate the impacting object, while the highly nonlinear structural response is analysed with the finite element method. Both methods are coupled within an open-source multi-physics research code to exchange data and simulate the interaction. The successful practical application of the coupling algorithm is demonstrated with this work, as the numerical results show good agreement with the experimental results. Within this paper the main focus is the appropriate modelling of the impacting objects, which heavily influences the simulation results, while a simplified structural model allows a correct assessment of the global deformation behaviour and reaction forces.


Introduction
Natural disasters, such as rockfall events, pose a serious threat to people, especially in mountainous areas. Several works, such as [1][2][3][4], already describe how to simulate these events and discuss the appropriate modelling of protection structures. While [5,6] use the coupling of the DEM and the FEM to analyse the interaction of rocks and the ground surface, [7] discusses different design strategies for flexible rockfall barriers.
Recently, the work of [8,9] explained how to realize a multi-physics analysis of rock impact in an open-source code. The code that is used within this work is KRATOS [10][11][12][13]. Since this is an open-source code, it allows for further development by independent research groups and is publicly accessible. The development and use of open-source codes allows the reliable assessment of life threatening natural events for everyone and is of upmost importance for the safety of the society.
To model the impact scenario two different numerical methods are combined in a coupled simulation. The finite element method (FEM) is used to analyse the structural response to the impact, while the discrete element method (DEM) models and simulates the impacting objects. While the DEM represents a particle method [14][15][16], the FEM is concerned with the analysis of continuous, deformable bodies.
In contrast to continuum-based particle methods such as the material point method (MPM) [17,18] or the particle finite element method (PFEM) [19,20], the DEM describes single discrete particles, instead of a continuum, which is composed of particles. In the case of debris flows and avalanches, a continuum-based particle method such as MPM or PFEM would prove to be more useful. In the case of the simulation of rockfalls, the modelling of the impacting objects as discrete objects seems to make the most sense.
The DEM models the dynamics of spherical particles taking into account external forces such as gravity and contact with other particles or objects [14][15][16]. It can be used for many different particle shapes such as rectangles, cones, spheres and more; however, in this publication only spherical particles are considered. To describe more complex shapes, a set of spheres is connected within clusters. A cluster consists of several particles, which are used for contact detection and force evaluation. The mass and centre of gravity are described within the cluster shape and independent of the masses of the particles (for more information see [14,15]).
Current state-of-the-art publications, such as [1,3], use the FEM to model the geometry of the test objects. By using geometric objects such as surfaces and lines, the use of FEM in this case leads to complex and costly contact algorithms [16]. The use of spherical particles greatly simplifies contact detection and thus appears to be the more appropriate modelling method. Additionally, the general restriction of simple FEM geometries is alleviated since particle clusters can represent arbitrarily shaped objects.
The objective of this work is to show the applicability of the aforementioned coupling method and discuss the appropriate modelling of impacting objects. In contrast to standard finite element discretization of the impacting objects as provided in previous works, e.g. [1][2][3][4], the DEM in KRATOS uses clusters of single spheres to represent any desired, arbitrary shape. This allows a simple contact search but also raises the question how these clusters have to be defined to achieve adequate results. Since this work focuses on the appropriate modelling of the impacting objects, the application of the coupling algorithm and not on the exact analysis of detailed structural components, the results of a test resulting in an elastic structural response are used. Although large rockfall events often result in much higher impact energies, often generating plastic and possibly destructive structural responses, this behaviour is not considered for demonstrative simplicity.
From a formal point of view, the structure of the paper is as follows: -Section 2 briefly describes the underlying coupling theory and the respective components in this multi-physics simulation (see [8,9] for detailed information). -Section 3 describes the experiment which will be used for validation and calibration. -Section 4 discusses the modelling of the numerical simulation used to replicate the experiment. Both the modelling of the structure itself and the modelling of the impacting object are discussed. -Section 5 presents the numerical results and explains each of the results in detail. -Section 6 finalizes this work with a conclusion and an outlook for future research.

DEM-FEM coupling
The general coupling of the two participating methods (DEM, FEM) was discussed in detail in recent publications such as [5][6][7]15,21]. In particular, [5][6][7] describe the possibility to apply the aforementioned coupling to analyse rockfall events modelling the interaction of the rock with the ground.
Recently [8,9] used the FEM to calculate the highly nonlinear structural response of protection nets to the extreme load cases of rockfall events (simulated by the DEM), allowing for a deeper understanding of the interaction of impacting objects with flexible tension structures. Nevertheless, for the sake of completeness, a short introduction is given in the proceeding subsections expressing tensors with bold letters and scalar values with regular letters. A list of appropriate, more detailed publications on the respective topics is attached to the specific subsections.

DEM
As a particle method, the DEM is used to model single discrete particles. The movement of discrete particles can be effectively simulated, making DEM an advantageous method for modelling impacting objects. It offers the advantage that rigid objects of any shape can also be modelled as particle clusters, which greatly reduces the computational effort. This feature is used in this work to model the impacting objects.
Various contact laws may be formulated which can be applied to the case at hand. The contact of particles to arbitrary geometric objects such as points, lines and surfaces is always modelled in the present implementation as the contact of the geometric objects to a sphere. This sphere has a radius R i and a geometrical centre C i .
For the present work three variants are of particular interest (described in [21]), where d e is used to express the smallest distance to the line element in contact. The geometric contact conditions are stated in Table 1.
Particle i and line ( Fig. 1) The corresponding contact forces are calculated with the contact law [21], which is suitable for the respective application, and thus, the equation of motion is integrated based on Newton's second law. In this study a Hertz-Mindlin springdashpot contact model (abbreviated as HM+D in [15]) is used. A sketch of the rheological contact model is shown in Fig. 2. The respective normal, tangent spring stiffness k n , k t as well as the normal, and tangent damping coefficients c n , c t are used to calculate the contact forces [9,15,[21][22][23]. To restrict the maximum tangential force to the Coulomb's friction limit the coefficient of friction μ is applied.
A detailed description of the contact algorithms, force evaluations, time integration, and other DEM-related topics would lie beyond the scope of this work. For a more detailed introduction to DEM see [14][15][16][21][22][23][24].

FEM
The FEM is used in the context of this work to discretize the protective structure. As described in [25], the virtual work is set up with the involvement of the d'Alembert forces and the external virtual work δW ext , As depicted in Eq. 1 the second Piola-Kirchhoff stress σ P K 2 , the work equivalent Green-Lagrange strain ε GL , the density ρ as well as the degrees of freedom u (here the displacements) and their second time derivativeü = ∂ 2 u ∂t 2 are integrated over the reference domain Ω 0 .
By evaluating Eq. 1 the following equilibrium equation results, with the mass matrix M, the damping matrix C, the stiffness matrix K, and the external forces (the DEM contact forces in this context) f ext . Because of the very short time interval investigated in the field of rockfall simulation [4], explicit time integration is suitable for solving this system of equations, since small time steps are required to resolve the impact process with sufficient accuracy. This method does not require the complex solution of sometimes very large systems of equations. A diagonal mass matrix is used (point mass assumption) and Eq. 2 is solved with the help of Newton's second law, using the internal forces f int and the discrete nodal velocitieṡ u.

Coupling procedure
The coupling procedure is described in detail in [8,9]. The basic idea is the exchange of data between two stand-alone analyses. Figure 3 describes this idea. The respective steps can be summarized in the following enumeration (numbers relate to Fig. 4). As known from other coupled multi-physics problems, such as fluid-structure interaction (FSI) [28], the direct explicit transfer of the interface data (forces, velocities, displacements) can lead to divergence problems in the staggered simulation. This problem is caused by large contact forces

Fig. 4
Strong coupling communication diagram [9] due to differences in velocities, acceleration, and highly different masses on both sides. In contrast to the weak coupling approach (direct exchange of data and subsequent proceeding in time [9]), the strong coupling approach, which is presented in [9] and depicted in Fig. 4, adds an additional iteration loop in each time step, which solves for the equilibrium between both numerical physics. This requires a Gauss-Seidel loop between DEM and FEM, which might need to be solved multiple times within one time step.
A more profound discussion of the aforementioned coupling scheme can be found in [8,9].

Experiment
The experiment set-up is described below and depicted in Fig. 5. In the course of this work, only the required comparative data will be published, while the remaining data are not publicly accessible as property of Geobrugg. A standardized (SAEFL) concrete block was used, with a total weight of 180 kg (the concrete block weighs in itself 175 kg, the attached wire rope strap weighs 5 kg bringing the total weight of the block to 180 kg) and an edge length of 0.41 m impacts a 3.9 × 3.9 m 2 DELTAX ® G80/2 [29] (see Fig. 6) net which is spanned into a CTI frame with rigid boundaries (see [2] for a detailed description). The rigid boundaries are in this case 5/8" shackles connecting the mesh to the steel frame.
The block is dropped from a height of 2.0 m which results in an impact velocity of v = √ 2.0 · 9.81 · 2.0 = 6.261484 m s . The same test is repeated twice, whereas (due to different pre-   [29] stresses) the initial sag, as a result of the dead load, varies from 0.05 to 0.10 m (labelled with exp_1, respectively, exp_2 hereafter). The tests are laid out to produce a rebound of the block without failure of any mesh wire.
Three main parameters are obtained from these tests by the means of a high-speed camera and an internal acceleration measuring device. The high-speed camera records at a resolution of 1280 × 1024 and a frame rate of 500 fps. The acceleration sensor is composed of a built-in triaxial accelerometer with a range up to 500 g and a sampling rate of 20,000 Hz and is incorporated in the test block. Deflection, i.e. the vertical displacement or sag of the mesh after impact, is measured directly through the analysis of the videos obtained from the high-speed camera. The block's velocity is calculated by following the block's trajectory over time in the high-speed video, based on the principle dx/dt. The velocity before impact is compared between the computed one, as explained above, and the result from the video analysis, to ensure the plausibility of the video analysis. The analysis of the accelerations obtained by the accelerometer yields the force evolution over time.
The repeatability of the tests is deemed to be given. Five tests were carried out in total. The test block and the drop

Structural modeling
Due to the set-up of the experiment and the fine mesh of the used protection netting (see Fig. 6) a simplification of the net to a closed, homogeneous surface suggests itself. This assumption simplifies the structural modelling and allows the use of two-dimensional finite elements. Publications such as [1,3] suggest a shell element to be used. While the advantages are clarified in the publication cited, this would mean an overhead for the underlying experiment. Due to its simplified set-up a membrane element [25] is used in this work. Carrying only in-plane stresses and omitting rotational nodal degrees of freedom an initially plane geometry possesses zero stiffness. As a remedy a minor pre-stress is applied to the structural model (see Table 2) to provide a non-singular stiffness matrix at the beginning of the simulation. The Young's Modulus is tuned to the given value in Table 2 in order to achieve an initial sag due to dead load of ≈ 0.05 m (see Fig. 7), as described in the experimental set-up.
With respect to the technical data sheets available on [29] the thickness of the membrane element is set to 0.008 m, whereas the density of the membrane is derived from the provided weight of 0.65 kg m 2 DELTAX ® mesh standard roll:  Fig. 7 Initial static analysis to obtain equilibrium position for dead load (plot scaled by factor × 10) As described in the experimental set-up in Sect. 3 the two tests that are investigated in this work show a rebound of the impacting sphere and no damage. Due to the observed structural response a simple elastic material model [25] and a Poisson's ratio of ν = 0 is applied in this work. In the course of this work it is shown that structural models with the right simplifications allow a correct assessment of the global behaviour. If a detailed analysis of the deformation behaviour and failure behaviour of individual structural elements is to be carried out, more attention must be paid to structural modelling. Publications such as [3,30] a contain a very detailed description of the modelling and the state of the art.

Impacting object
In contrast to preceding works, such as [1,2], this work proposes to model arbitrary objects with the help of discrete spherical element clusters. The advantage over standard finite element discretized objects is the simplified contact detection between arbitrary boundary objects and single spheres contained in the cluster [21].
As depicted in Fig. 8, the standardized experiment object is modelled with seven levels of refinement. Ranging from a representation as a single sphere to a detailed geometrical description with up to 22,232 spheres. Special algorithms are needed to create such refined cluster files. This work uses the sphere-tree algorithm described in [31,32] as it is available in an online toolkit [33].
The density of the respective cluster is fitted to obtain a total cluster mass of 180 kg (see Sect. 3). Important DEM parameters, such as the coefficient of restitution ε [22] and the Young's Modulus of the particle E p , are varied with respect to Table 3 with their influence studied. The range in which E p is selected is based on experience gained from previous simulations.

Validation
The experiment described in Sect. 3 produced results which are used to validate the proposed coupling algorithm [8,9] and to investigate the influence of discrete element cluster  Each of the aforementioned clusters, described in Table 3, is used in the numerical impact simulations. If not mentioned otherwise, a time step value of d t = 10 −5 s is used. The results with the label "small dt" follow from a simulation with d t = 10 −6 s. Although the simulated system set-up is more similar to exp_1, the experimental results for exp_2 are added to the result plots as it represents an experiment with the same impact velocity and the same of mass of the impacting object. The difference in the initial sag is a results of the test set-up and can hardly be controlled; thus, it allows error zones to be added in the result plots to show the variability of the experimental results.

Displacement
The simulated displacement values of the impacting object are compared to the experiment results. Figure 9 shows that the simulation results lie between the two experiment results.
It is noticeable that the refinement of the cluster shows a large influence as soon as the number of spheres is larger than one (c 1 ). To explain this behaviour the simulation is visualized at the time of maximum displacement in Figs. 11, 12, 13, 14, 15, 16 and 17 (see Fig. 10 for comparison to the experiment). Due to the fact, that the contact forces in DEM are calculated with the help of spring-dashpot models [15,22,23] the single sphere with one large radius experiences an indentation which is larger than it is for the clusters with more (smaller spheres). This has additional effects too, which become clear in Figs. 18 and 19. The single particle has a larger indention (Fig. 11) and thus decelerates slower (Fig. 18). This finally leads to lower contact/reaction forces and a longer duration of load application, as visualized in Fig. 19. The same behaviour is observed when decreasing the time step (see Fig. 9, label "c1 small dt") and thus can be concluded to be a problem of the very rough representation of the original geometry as a single sphere.
Regarding the remaining refinement levels, a good agreement with the exp_1 is shown, as the initial sag of the numerical model due to dead load matches with the initial sag of the first experiment set-up exp_1 of 0.05 [m], while the displacements of all simulation results are within the error zone.

Velocity
Similar to the observations in Sect. 5.1 all cluster refinement levels again show good agreement with the experiment results in Fig. 18. The exception is again the single sphere which decelerates more slowly due to its larger indention (Fig. 11).

Forces
Lastly, the forces obtained from the numerical simulation are compared to the experiment results. Figure 19 shows the reaction forces in the structure as a sum over all boundary nodes. A good agreement for all cluster refinements >c_1 can be observed for both the force value as well as the contact time. Similar to the observations in the previous sections the roughest object representation performs poorly. The maximum force is below the experimentally obtained value and the time in which the object is in contact with the structure is too large. c_1, as the roughest object representation, is a single sphere, poorly representing the correct object shape. With regard to the forces, it can also be said that a more precise modelling of the exact object geometry plays a major role, even if the impact, as in this case, takes place without rotations.
One more interesting behaviour that can be observed in the simulations is the difference between reaction forces F R (structural response, FEM) and contact forces F C (cluster  interaction with boundary, DEM). Figure 20 illustrates the difference for a varying time step size, which is shown in Table 4. For simplicity, cluster c1 is used. The theoretical difference should be described by only the dead load of the protection net. Nevertheless strongly oscillating contact forces are observed for rather larger time steps, while the reaction forces remained smooth.
As visualized in Fig. 20, the contact forces F C for the largest time step (s1) strongly oscillate, while the respective reaction forces F R develop smoothly and exactly in the middle of the oscillations. This behaviour can also be observed for the next smaller time step (s2) as soon as the maximum force is reached. For the smallest time step (s3) the described behaviour is almost negligible.
Considerably large contact forces in combination with large time steps that are too big to resolve the contact interaction can lead to particles losing contact in one time step only to regain contact again in the next. As a result of this fretting, oscillating contact forces F C are generated. Since the contact points are located locally in the centre of the structure and are at some distance from the edge nodes at which the reaction forces are calculated, the influence of the initial oscillations

Results
With respect to the relative computational time for each simulation, as shown in Fig. 21, the effect of additional particles is negligible under ≈ 10 3 spheres. Additionally, considering the results in Sect. 5, the cluster version c5 is recommended for usage. It represents the best compromise between accuracy and computing effort, while also properly representing number of spheres in cluster [-] Fig. 21 Comparison of relative computational time the geometry (see Fig. 8). It is also expected that the fine representation of the impacting object influences the simulation results favourably if the rotation of the object plays a crucial role (see also [8], where a simplified Attenuator is simulated by using an arbitrary shaped rock object). The computational time for the simulation of the cluster c5 is discussed as an example. A total time of 1108s ≈ 18.5 min with an increased time step of d t = 2 · 10 −5 s for a total simulation time of 0.3 s is required. With respect to Fig. 9, the cluster leaves the reference axis of origin at ≈ 0.25 s; thus, the simulation time could be further reduced by ≈ 17%. The CPU system settings for this simulation are an Intel(R) Xeon(R) CPU E5-2623 v4 @ 2.60 GHz.

Conclusion
This contribution demonstrates that using the DEM for modelling an impacting object as well as simulating the interaction with the FEM structure is expedient. It allows modelling of any desired shape without requiring complex contact models between surfaces, lines and vertices. The strategy of using clusters of single spheres simplifies these contact models to a simple sphere contact detection. As shown in this work, this feature can be used to simulate standardized test objects, as illustrated in Fig. 8, as well as arbitrary shapes (see Fig. 22).
In cases of real rockfall scenarios it can be beneficial to simulate different, varying rock shapes to also capture the influence of sharp edges, and progressive fragmentation which result in a reduction of kinetic energy [6]. Additionally smaller sizes of impacting objects can cause a bullet effect [4]. The use of discrete element clusters in combination with cohesive forces, as described by [15], allows the cluster to break, reaching limit stresses, and ultimately allows modelling of cluster fragmentation.
Additionally a detailed discussion of the comparison between simulation results and test results highlights that a coarse description of the impacting object renders poor results (Figs. 9, 18, 19), while a reasonable representation of the impacting object geometry can be achieved without significantly increasing computational time (Fig. 21). The simplification of the structure to a closed surface (see Sect. 4.1) and the appropriate modelling of the impacting object (see Sect. 4.2) are confirmed to strike a good compromise between computational time and the approximation of the physically correct structural response (see Figs. 9,18,19). The relatively simplified modelling of the structure was shown to predict reasonably accurate global deformation behaviour and reaction forces. In combination with the highperformance modelling of the impacting objects as spherical clusters, this approach has advantages over current publications that use very complex structural models and thus Fig. 22 Arbitrary rock shape modelled with discrete element cluster require a lot of computing time. The contact forces and reaction forces of the structure were also compared. It was shown that relatively large time steps can lead to oscillating contact forces, which can be counteracted with smaller time steps restricting its influence on local parts of the structure. Although the contact forces can oscillate strongly, the reaction forces are relatively smooth. In addition to the elastic tests considered in this paper, experiments that lead to plastic deformations or even damage to the mesh were carried out too. It is expected that the application of plasticity constitution laws and damage models in the FEM part will also allow the reproduction of the results of the other experiments.
The accurate assessment of life threatening natural events (such as rock slides, mud flow, and strong wind events) is of great value to society, with the open-source KRATOS [10,12,13] multi-physics software offering advanced analysis capabilities free of charge for all. The code is publicly accessible via a GitHub repository [34] with an installation guideline provided. KRATOS is designed in C++ and includes a Python interface to facilitate the advanced development and simulation. With easy access and a community where participation is highly appreciated, there is a great potential to further advance these simulation technologies. Such open-source projects are of upmost importance for the society, especially in times of extreme weather and environmental events.
Author contributions All authors prepared the manuscript. All authors read and approved the final manuscript.
Funding Not applicable.

Availability of data and materials
The comparative data in this work are the results of several experiments conducted by Geobrugg. The experiments were carried out in 2018 in Walenstadt, Switzerland, according to the Swiss guideline (SAEFL). The data are the property of Geobrugg and are partly published in the course of this work; the remaining data are not publicly available.

Compliance with ethical standards
Competing interests The authors declare that they have no competing interests.

Code availability
The software used is [10,12,13]. The current developers' version is available at [34] as a GitHub repository.
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://creativecomm ons.org/licenses/by/4.0/.