Advanced Modeling and Simulation of Rockfall Attenuator Barriers Via Partitioned DEM-FEM Coupling

Attenuator barriers, in contrast to conventional safety nets, tend to smoothly guide impacting rocks instead of absorbing large amounts of strain energy arresting them. It has been shown that the rock’s rotation plays an important role in the bearing capacity of these systems. Although experimental tests have to be conducted to gain a detailed insight into the behavior of both the structures and the rock itself, these tests are usually costly, time-consuming, and offer limited generalizability to other structure/environment combinations. Thus, in order to support the engineer’s design decision, reinforce test results and confidently predict barrier performance beyond experimental configurations this work describes an appropriate numerical modeling and simulation method of this coupled problem. For this purpose, the Discrete Element Method (DEM) and the Finite Element Method (FEM) are coupled in an open-source multi-physics code. In order to flexibly model rocks of any shape, sphere clusters are used which employ simple and efficient contact algorithms despite arbitrarily complicated shapes. A general summary of the FEM formulation is presented as well as detailed derivations of finite elements particularly pertinent to rockfall simulations. The presented modeling and coupling method is validated against experimental testing conducted by the company Geobrugg. Good agreement is achieved between the simulated and experimental results, demonstrating the successful practical application of the proposed method.


INTRODUCTION
Rock impact is an exceptional load case due to many factors. The shape, the density and the motion process itself of the impacting object can be arbitrary and difficult to predict. The additional dynamics and arbitrariness of the impact point make it impossible to determine a generally valid load assumption formula. This raises the question of how to design protective structures for such load cases. During rock impact events, protective structures are generally engineered to deform greatly, which reduces peak system decelerations and safely dissipates large amounts of energy.
In order to allow the protective structure to undergo large deformations, a wide variety of details (such as ring elements and braking elements in Figure 1B) are installed, the direct finite element modeling of which represents a demanding task.
Rockfall protective structures may be divided into two different categories: active and passive measures. While nearsurface nets (see Figure 1A) for slope stabilization represent a so called active measure and are easier to dimension, passive protective measures are installed as soon as the cause of the rockfall cannot be prevented. Passive protection structures are the focus of the present work and can in turn be divided into selfcleaning and non-self-cleaning (see Figures 1B,C). In the course of this work an experiment of an Attenuator barrier protective structure (illustrated in Figure 1D) is numerically replicated. The principle of an Attenuator is to intercept rockfall trajectories and attenuate their kinetic energy, while guiding the rockfall to a designated catchment area. These systems have been employed for several decades as described in Muhunthan et al. (2005) and were the subject of many field tests (Arndt et al., 2009;Glover and Ammann, 2016;Wyllie et al., 2016). An Attenuator barrier rockfall experiment program conducted by the company Geobrugg in 2017 is the primary experimental focus in this work.
Experiments confirm that the rotation of the impacting object has a substantial influence on the performance of Attenuator barriers. This behavior is further investigated with the help of numerical simulations in the present work. Throughout the experimental program, an accelerometer was installed in the center of the test objects and recorded data such as accelerations and angular velocities. Since the exact orientation of the rock at impact cannot be determined precisely, these data are only compared in terms of their general trends instead of exact values (see section 5.2). In order to better understand the interaction between impacting objects and highly flexible protective structures and to better answer practice-oriented questions, numerical methods are used. These allow a detailed insight into the structure, its utilization and have the potential to show optimization possibilities. By employing coupled simulations, different numerical methods, each with their individual strengths and domain of application, may be brought together to model and analyze the interaction of different physical problems.
The present work couples the Discrete Element Method (DEM) with the Finite Element Method (FEM) to simulate the interaction between impacting objects and highly flexible protective structures. The DEM is used to simulate the impacting object while the FEM is used to calculate the appropriate dynamic structural response. The advantages of this choice of simulation methods are explained in the following sections. To demonstrate the accuracy and applicability of the aforementioned coupling strategy the 2017 Geobrugg experimental program is modeled and simulated. The resulting data show good agreement with the experimental data demonstrating the efficacy of the numerical method, as well as confirming the effectiveness of the Attenuator barrier design. Two different experiments are simulated to guarantee the validity and applicability of the method presented here.
The discrete element method is a discrete particle method with unique strengths in analyzing the movement and interaction of individual particles. The method was first proposed by Cundall and Strack (1979) and has been increasingly adopted since. Matuttis and Chen (2014) describe in detail the underlying theories of DEM while many publications deal with the precise analysis of contact formulations. Special mention should be made of Shäfer et al. (1996) who has done some fundamental work in this topic. Building on this, Schwager and Pöschel (2007); Cummins et al. (2012); Thornton et al. (2012) have described improved formulations and contact models, contributing to the current state of the art. A detailed description of contact search and calculation of contact between discrete spherical particles and boundary conditions discretized by finite elements was provided by Santasusana (2016), Santasusana et al. (2016) which lay the foundation for the present work. Gao and Feng (2019) extended this idea by exchanging the FEM with the isogeometric method, however, this is not expected to confer any advantages for this paper's scope. Irazábal et al. (2019) conducted a detailed investigation of suitable time integrators, especially with regard to the integration of the rotational motion, which is of great importance in the present case. Specifically related to the case of rockfall Chau et al. (2002) studied the appropriate coefficient of restitution for rockfall events. The use of simple spherical particle clustered together to describe arbitrarily shaped objects in DEM was formulated by Madhusudhan et al. (2009) and subsequently employed by Sautter et al. (2021) for the analysis of rotationfree rockfall events. This flexible modeling strategy, employed in this paper, allows the use of simple contact algorithms without sacrificing the approximate modeling of irregular rock geometries. As already mentioned, experiments confirm that the rotation of the impacting objects has a large influence on the overall behavior of the impact, which is why the shape of the rock has to be taken into account. Disc shaped one for example have normally much higher rotational energy compared to cube shaped blocks.
The FEM is used to numerically solve partial differential equations to analyze continuous, deformable bodies. For the interested reader, notable FEM references include Basar and Weichert (2000), Belytschko et al. (2000), Holzapfel (2010). A general FEM summary follows in subsection 2.2 that also includes a detailed derivation of the finite element formulations used in this work.
The coupling of the finite element and discrete element numerical methods builds on the work of Santasusana (2016), Santasusana et al. (2016) and was continued in Sautter et al. (2020). The implementation is in an open source multi-phyiscs code called KRATOS based on the work of Dadvand et al. (2010), Dadvand et al. (2013);Ferrándiz et al. (2020). The importance of an open source solution for the analysis and assessment of natural events such as rockfall [studied in detail by Volkwein (2004), Volkwein (2005); Zhao et al. (2020b)], strong wind [fluidstructure interaction described by Wüchner (2006)] and debris flows [discussed in detail by Wendeler (2008), von Bötticher (2012), Zhao et al. (2020a)] cannot be overemphasized. In combination with the detailed description of element formulations, contact algorithms, and coupling methods this work offers the possibility to carry out independent analyses of any rockfall events and the evaluation of suitable protective structures.
In principle, it is conceivable to employ other particle methods for the analysis and simulation of the impacting objects. In contrast to DEM, the Material Point Method (MPM) and the Particle Finite Element Method (PFEM) are continuum-based particle methods. Due to their continuum-based approach, both methods possess particular strengths compared to DEM, for example, in the modeling of debris flows or avalanches. Due to the properties of rockfall, DEM is preferable in this case [see also Lisjak and Grasselli (2010), . Further information on MPM can be found in Zhang et al. (2016), Chandra et al. (2019), Larese et al. (2020), Wilson et al. (2020), while Salazar et al. (2015), Larese (2016) describe PFEM.
A comprehensive review and prospective on the subject of the simulation of rockfall protection systems is given in Volkwein (2004), Volkwein (2005) and represents the basis for the current state of the art. The work lead to Escallon et al. (2015) describing highly detailed modeling of flexible protection structures and their components. The exact modeling of the individual components often turns out to be overhead, especially when the global structural behavior is of interest. In order to reduce the complexity of the simulations while accurately analyzing protection structure components at selected local areas, a hybrid model is described in Tahmasbi et al. (2019). The idea is to use a simple element formulation in the boundary regions of the protective structure while local locations of interest are modeled with the formulations from Escallon et al. (2015). Other works by Sasiharan et al. (2006), Dhakhal et al. (2011), Mentani et al. (2018 among others, also employ the concept of representing the complex protection structure with a simpler surface description composed of shell and isotropic membrane elements (Tahmasbi et al., 2019). In the present work, we will take up the idea of membrane elements, and incorporate the anisotropic material law of Münsch and Reinhardt (1995) to represent the different behavior in the two main directions that can be observed in experiments.
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. A summary and description of the finite element formulations used within this work is added in subsection 2.2. • Section 3 discusses the 2017 Geobrugg experiment used for validation and calibration. • Section 4 presents the modeling of the numerical simulation used to replicate the experiment. Both the modeling of the structure itself (in subsection 4.1) and the modeling of the impacting object (in subsection 4.2) are covered. • Section 5 illustrates the numerical examples and their results across displacements (in subsection 5.1), translational velocities (in subsection 5.3), and angular velocities (in subsection 5.2). • Section 6 finalizes this work with a conclusion and an outlook for future research.

DEM-FEM COUPLING
To realize the coupled simulation, the DEM and the FEM are coupled in a staggered framework, summarized in subsection 2.3 [refer to Sautter et al. (2020) for full details], which allows the use of any two standalone solvers for the DEM and the FEM components. In principle, many different methods can be used to analyze the problem at hand. In order to find the combination of suitable methods, one has to think about the general characteristics of the problem at hand. Rockfall is a discrete event, so one must be able to model discrete objects. This can be solved with the FEM, as shown by Volkwein (2004), Volkwein (2005), Sasiharan et al. (2006), Dhakhal et al. (2011), Escallon et al. (2015; Mentani et al. (2018), Tahmasbi et al. (2019). The surface of the discrete objects is modeled with the help of points, lines, and surfaces. If one now wants to find the contact to a structure that is also modeled with points, lines, and surfaces, very complex contact algorithms are necessary. The contact between surface/surface, surface/line, and line/line is sometimes very demanding and computationally intensive. If, on the other hand, a particle method is used that represents/approximates the discrete objects as simple spheres, the contact check is extremely simplified, see the work by Santasusana (2016), Santasusana et al. (2016). Contact algorithms for the combination sphere/surface, sphere/ line, and sphere/vertex are very effective and fast to perform. By using clusters of spherical particles, arbitrarily shaped objects can be modeled despite the efficient contact algorithms for spheres (Sautter et al., 2021) describes the successful application of these clusters already for the simulation of other types of rockfall protection systems]. To decide on the appropriate particle method, it is again helpful to look at the discrete characteristics of the rockfall event. Other particle methods such as MPM [first proposed by Sulsky et al. (1994)], PFEM (Cremonesi et al., 2020), or Smoothed-Particle Hydrodynamics(SPH) [first proposed by Gingold and Monaghan (1977), Lucy (1977)] are particularly suitable for the simulation of non-discrete natural events such as debris flows, avalanches, or shallow landslides. For this reason, DEM is chosen to model and simulate the impacting objects (in this case rocks). To model the structure the FEM is used, as it has proven to be very effective and well established. It allows to model special structural elements and to effectively assess structural stress and deformation states. With the help of the coupling environment in KRATOS [see Dadvand et al., (2010), Dadvand et al., (2013), Ferrándiz et al. (2020)], which we co-developed, these methods can be efficiently combined with each other, providing a suitable simulation environment for the problem at hand. Subsection 2.1 briefly describes the DEM and introduces key aspects to give a general overview of the topic. The FEM is briefly summarized in subsection 2.2 which includes detailed descriptions of the element formulations employed in this study. Throughout this work, bold symbols indicate tensors while italics denote scalars.

DEM
In contrast to continuum-based particle methods like the MPM (Zhang et al., 2016;Chandra et al., 2019;Larese et al. 2020;Wilson et al. 2020]) or the PFEM, the DEM is a numerical method to simulate the motion and interaction of discrete particles. Since Cundall and Strack (1979) introduced the first considerations regarding DEM, the method has been increasingly adopted for the analysis of discrete particles both in research and industry. The motion and interaction between particles and also between particles and finite element discretized boundary conditions can be simulated very effectively, see Santasusana (2016), Santasusana et al. (2016). In the present work only spherical particles are used which reduces the effort of contact search with arbitrary geometries enormously and thus leads to a very efficient simulation method. Although spherical particles are unable to model complex geometries by themselves, this work combines spherical particles together into clusters to model complex rock geometries while still using the simple contact search of the spherical particles. Madhusudhan et al. (2009) has studied these in detail and Sautter et al. (2021) has already used clusters of spheres for simple rockfall events. Chapter 4.2 describes in detail how this strategy is used in the simulation.
The search for contact with spherical particles has been studied in detail and is exhaustively discussed in Santasusana (2016), Santasusana et al. (2016). Each particle has a radius R i and a geometrical center C i . Contact is found when the distance between C i and another object (particle, line, area or vertex) is smaller than the respective radius R i .
As soon as contact is detected, contact forces are calculated. There are a number of different contact laws available, and the user has to select the most appropriate for the current case at hand. In this work a Hertz-Mindlin spring-dashpot contact model was used (abbreviated with HM + D). The corresponding rheological model is shown in Figures 2A,B. To calculate the contact force, the normal and tangential spring stiffness k n and k t , the normal and tangential damping coefficients c n and c t , and the friction coefficient μ must be determined. μ is used to limit the tangential force to Coulomb's friction limit.
Finally the equation of motion according to Newton's second law is integrated in time and one advances to the next time step.
More detailed descriptions of the DEM can be found in Santasusana (2016), Santasusana et al. (2016). Matuttis and Chen (2014) provide an overview of the whole method and also discuss the use of non-spherical particles. A detailed discussion of appropriate time integrators may be found in Irazábal et al. (2019) while Schwager and Pöschel (2007), Cummins et al. (2012), Thornton et al. (2012), Shäfer et al. (1996) discuss the modeling of non-elastic collisions and the appropriate choice of the coefficient of restitution.

FEM
The FEM is used in this work to discretize the protective structure. As described in Belytschko et al. (2000), Holzapfel (2010), the virtual work δW is set up with the involvement of the d'Alembert forces, the external virtual work δW ext , and the internal virtual work δW int , As depicted in Eq. 1 the second Piola-Kirchhoff stress S, its energy-conjugate Green-Lagrange strain ε, the density ρ as well as June 2021 | Volume 7 | Article 659382 the degrees of freedom u (here the displacements) and their second time derivative € u z 2 u zt 2 are integrated over the reference area Ω 0 . This equilibrium is solved numerically with a Newton's type iterative solution technique operating on linearized virtual work components as discussed in Belytschko et al. (2000). Internal element forces F int,r are expressed with respect to each degree of freedom r in the system, according to, Additionally, the tangent stiffness matrix K, is often needed for Newton's type solution schemes [see also Belytschko et al. (2000), Holzapfel (2010)].
In the subsequent sections lower case symbols represent quantities in the current configuration (x) while capital letters describe quantities in the reference configuration (X).
In view of complex protective structures, such as Attenuator barriers, the structure must be modeled accurately with attention to detail. This work introduces the application of two unique finite element formulations to be used in the simulation of such structures: the sliding cable element and the plane-stress membrane element.

Sliding Cable Element
Compliant structures, such as Attenuator barriers rely on the ability of large deformations to reduce peak decelerations and withstand impacting objects. One feature which is used to realize large deformations is the "loose" connection of the net to upper ropes as illustrated in Figure 3. The net is connected orthogonal to the dark blue support rope but is allowed to slide along the rope, impeded only by friction.
To properly model this behavior in a FEM simulation a sliding cable element formulation [Volkwein (2004), Boulaud and Douthe (2017)] is employed, which also enables to integrate a deflection roller into the simulation (left side of Figure 3A).
Instead of single linear cable elements connecting the nodes n i and n i+1 (see Figure 3B) the element formulation used in this work considers all discrete nodes along a predefined path in one single element. This is realized by computing the total strain of the whole cable element with respect to the total length l instead of the single lengths l i . According to Volkwein (2004), Boulaud and Douthe (2017) the total length l in the current configuration and the total length L in the reference configuration are calculated as the sum of the respective line segments, As described in Belytschko et al. (2000), Holzapfel (2010) these quantities are used to express the 1-dimensional Green-Lagrange strain ε, and the 1-dimensional 2 nd Piola-Kirchhoff stress S, with the help of the Young's Modulus E. Including the prestress S pre the virtual internal work δW int as stated in subsection 2.2 is expressed, where A represents the reference cross-section of the cable element. Following the derivation in Boulaud and Douthe (2017) an explicit formulation for δW int is developed, Equation 7 describes the virtual internal work of the sliding cable element with the help of the ratio between the total current length l and the total reference length L, including the vector of nodal virtual displacements δu and the direction of internal forces at each discrete node of the element assembled in vector T, Here n denotes the total number of nodes and the coordinates x, y, z describe the current configuration and the nodal coordinate distances, e.g., Δx i x i+1 − x i . Additionally, the friction is included and needs to be evaluated at each discrete node in the element by first calculating an equal and opposite force F i F int,i and subsequently using this force to calculate a friction force ΔN i μF i according to the friction coefficient μ. Volkwein (2004) discusses that this force must now be considered in addition to the normal elastic forces in the rope.
To validate the correct implementation of the sliding cable element, a simple gravity driven impact simulation is illustrated in Figures 4A,B, demonstrating the expected large sliding deformation of the net shackles.
Alternative solutions to enable the sliding of nodes along a cable element can be achieved using the penalty method [see Bauer et al. (2018)], Lagrange multipliers [see Holzapfel (2010)], or multi point constraints (MPC) as described in Abel and Shephard (1979). While the first two methods minimize the orthogonal distance between an arbitrary node and the axis of the sliding cable, the MPC approach directly couples degrees of freedom in such a way that the node travels along the desired path. Compared to these abstracted approaches, the sliding cable formulation, which inherently represents sliding nodes in its formulation, is the most robust solution to the problem presented and adopted in this work.

Plate in Membrane Action
As described in Sasiharan et al. (2006), Dhakhal et al. (2011), Sautter et al. (2021 membrane finite elements are suitable to approximate the behavior and response of rock-fall cable nets. With respect to the mesh geometry illustrated in Figure 6C, an anisotropic material law defined in Eq. 11 appropriately describes the net's direction-dependent behavior and is employed henceforth.
Derivation of the 2-dimensional plane-stress membrane element within 3-dimensional space commences by defining the general covariant base-vectors that span an arbitrary angled local coordinate system {G i } and {g i } expressed by, in which the standard FEM shape functions N i are employed (Belytschko et al., 2000). As the element is always a 2-dimensional plane within a 3-dimensional frame, the base-vectors G 3 and g 3 are constructed as a unit vector normal to that plane. Belytschko et al. (2000) notes that these arbitrary base-vectors are not necessarily orthonormal and are used to calculate the Green-Lagrange strain tensor ε, where G ij and g ij is the metric of the respective configuration and G i expresses the reference contra-variant base-vectors. In order to express the elastic stress-strain relation S C : ε in clear matrix notation the stress and strain are arranged in Voigt notation {•} yielding {S} C · {ε}. C describes a 4 th order material tensor, whileC is of 2 nd order.
The anisotropic elastic consistent linearized tangent modulus C from Münsch and Reinhardt (1995) dependent upon the Young's Moduli E x E y , Poisson's ratios ] xy , ] yx , and the shear modulus G is expressed as, Since the {G i } basis is not necessarily orthonormal, the strains naturally defined in that basis must be transformed to a local orthonormal basis {A i } {A i }. The construction of such a basis is achieved by normalizing G 1 and G 3 creating A 2 such that it is orthonormal to A 1 and A 3 , FIGURE 4 | Sliding cable elements on two opposing sides applying the tension field theory, described by Nakashino and Natori (2005), to exclude compression stiffness from the plate in membrane finite elements. (A) Reference configuration, fixed on four corner nodes. (B) Sliding cable element formulation allows for a sliding of the discrete element nodes.
Frontiers in Built Environment | www.frontiersin.org June 2021 | Volume 7 | Article 659382 The covariant strain coefficientsε ij in the local orthonormal basis {A i } transform by the following rule, given in Kiendl (2011),

Coupling Procedure
A partitioned coupling approach is employed to couple independent discrete element and finite element methods together in this work. Felippa et al. (2001) note that partitioned solution schemes possess the notable benefits of software modularity and also application of the most appropriate discretization and solution techniques to each individual component. Although this approach allows the use of arbitrary DEM and FEM codes, the user must consider the coupling of the two. Santasusana (2016) has already discussed this in detail and Sautter et al. (2020) has continued this coupling by introducing an additional Gauss-Seidel loop on the interface Γ level to ensure the convergence of the interface. The basic idea of the coupling method is the exchange of data between the two independent solvers. The solver of the DEM calculates contact forces F contact which are dependent on the displacements u and velocities _ u of the particles P and the discretized boundary conditions Ω D . These forces are mapped [see e.g., Tianyang (2016)] to the finite element mesh of the structure Ω S and by means of FEM the displacements and velocities are calculated which are finally mapped to the boundary conditions of the DEM solver. Multiple studies Wüchner (2006), Winterstein et al. (2018), Sautter et al. (2020) have confirmed that the simple exchange of this data can lead to divergence of the interface conditions, all of which are dependent on the time t. Figure 5 describes the use of additional iteration loops at each time step to achieve the convergent strong coupling as described by Sautter et al. (2020) (which also includes more details of the coupling). The steps in Figure 5 are summarized below:

2017 GEOBRUGG ROCKFALL EXPERIMENTS
In 2017 the Swiss company Geobrugg conducted an experimental program to confirm the effectiveness of a novel SPIDER mesh (Geobrugg, 2020). The following section is a brief summary of key test aspects relevant for numerical replication, refer to Hofmann et al. (2019) for full testing details.

Test Site
The test site was a near-vertical approximately 60 m high slope situated in a disused granite quarry in British Columbia, Canada. The Attenuator design investigated was the hanging net style with limited slope contact of the netting as per Figure 6A. The test set up consisted of three 8 m steel posts set 9 m apart and anchored to the slope with six retention ropes, two lateral support ropes, and one top rope penetrating through the post heads. The netting is attached to the top rope with shackles in each mesh opening as detailed in Figure 3A. Three principal zones of attenuation were investigated by varying the impact, transition, and collection zones of the high tensile steel wire nets. The experiments of 2017 are described in greater detail in Hofmann et al. (2019).

Site Survey and Rockfall Modeling
A terrain survey of the site was performed by Truline Survey in 2016 which was used as basis for rockfall modeling to adapt the Attenuator geometry between the testing series of 2016 and 2017. Particularly relevant to this paper, Sautter et al. (2020) demonstrates how such terrain models may be incorporated into simulations to further enhance accuracy and predictive confidence.

Impacting Object
The rocks used for testing were either granite blocks of various characteristic shapes (such as cubic, angular, and disc-like) or prefabricated concrete blocks with a special housing in the center for rock motion sensors (see Figure 6B). The blocks and rocks were then released by an excavator on top of the slope. The rocks bounced 3-4 times on the granite slope before impacting the netting.

Instrumentation
Experimental instrumentation included load-cells in all ropes, high-speed cameras filming the block's trajectory from different angles and rock motion sensors in the concrete blocks. This work will predominantly focus on the trajectory path, translational velocity and rotational velocity experimental data collected.

High Speed Cameras
Analysis of the high speed camera videos with a sampling rate of 500 fps facilitated quantification of the rockfall dynamics. By tracking the position of the block through time a translational velocity could be obtained and the tracking of every 90 + rotation of the block served as an indicative quantification of angular velocity. The cameras were located so as to have a frontal view on the mesh and a perpendicular view on the mesh. This allowed also to determine in which plane the block was moving.

Rock Motion Sensors
The rock motion sensors measured the accelerations and rotations of the block about its three axes for the duration of each test. A comparison of the angular velocity obtained from the video analysis and the rock motion sensor shows good agreement and therefore provides confidence in the translational velocity data estimated by video analysis.

Tests Used for Numerical Validation
The tests used to validate the modeling and simulation used a Geobrugg SPIDER S4/130 mesh [described in Geobrugg (2020) and visualized in Figure 6C]. Additional weights in the form of steel bars were shackled to the bottom of the mesh (see Figure 6A) to increase the inertia and vertically pre-tension the system and is incorporated in the simulation with additional point masses (discussed in subsection 4.1.1). The middle post was slightly bent after sustaining two rock impacts in previous tests, but the integrity and function of the system was not compromised, accordingly testing was continued, but it was chosen to incorporate it anyway in the modeling of the structure. The blocks used during the test, were a near-cubic 626 kg concrete block with an initial volume of approximately 0.75 × 0.75 × 0.75m 3 called T092 and a 278kg granite block with an initial volume of approximately 0.75 × 0.51 × 0.48m 3 called T089. Since the blocks were used for several tests and accumulated greater damage after each run (see Figure 6B for T092), the mass before and after each run was the characterizing measurement (instead of volume). Key experimental data pertinent to numerical replication including translational and rotational velocity at the time of netting impact and the rock's path and mass are summarized in subsection 4.2. Figure 7 is adapted from Geobrugg (2019) and depicts the respective participants in the simulation, each of which are discussed below.

Structural Modeling
Due to the small time steps necessary to resolve the impact, an explicit time integration scheme was selected (instead of an implicit scheme). The central-difference explicit scheme as described by Belytschko et al. (2000) was used to conduct the numerical time integration of the structural response with a time step of dt 5 · 10 − 5 s.

Point Masses
The additional weights at the bottom of the protection net, depicted in Figure 6, are modeled with single point masses (summing to an additional total mass of M add 358kg) equally distributed along the lower boundary of the membrane mesh (see Figure 7). Including these masses is critical to properly model both the gravitational forces and the additional dynamic inertia.

Sliding Cable Element
Modeling of the upper support rope, illustrated in Figure 3, demands a sophisticated element formulation, which is discussed in subsection 2.2.1. This formulation allows internal nodes to slide along the element, while internal forces are calculated based on the change of the total length, as discussed in Volkwein (2004); Boulaud and Douthe (2017). The following structural properties, based on construction plans in Geobrugg (2019), were applied in the simulation: Young's Modulus E 1.130 · 10 11 [N/m 2 ], cross area A 1.645 · 10 − 4 [m 2 ], density ρ 7.850 · 10 3 [kg/m 3 ], and the friction coefficient μ 0.25.

Standard Cable Element
Bracing cables, modeled with standard cable elements, are anchored in the rock and connected to strategically selected points on the supports of the protective structure. As presented in Belytschko et al. (2000), a 1-dimensional truss element formulation is applied and combined with an additional check for compression stresses in the element. If such stresses are detected the elemental stiffness contribution is temporarily eliminated from the global system of equations. Realistic structural properties were also obtained from construction plans provided in Geobrugg (2019) and listed in the following: Young's Modulus E 1.100 · 10 11 [N/m 2 ], cross area A 1.160 · 10 − 4 [m 2 ], and density ρ 7.850 · 10 3 [kg/m 3 ].

Posts
The posts, which predominantly carry vertical loads to support the protective structure, are simply supported in the rock. At the top they are connected to both the upper support rope as well as the bracing cables (seeFigure 3). In accordance with Belytschko et al. (2000), these posts may be suitably modeled as 1dimensional truss elements (neglecting any initial damaged bending), with the following data: Young's Modulus E 2.100 · 10 11 [N/m 2 ], cross area A 2.010 · 10 − 3 [m 2 ], and density ρ 7.850 · 10 3 [kg/m 3 ].

Protection Net
Due to the set-up of the experiment, similar to the set-up in Sautter et al. (2021), and the complicated geometry of the SPIDER net system (see Figure 6C), the net has been idealized as a closed homogeneous surface discretized with plane-stress membrane finite elements. Publications such as Mentani et al. (2018), Tahmasbi et al. (2019) suggest shell elements to be used, although they introduce additional complexity and computational expense compared to membrane elements. In the interest of simplicity, membrane elements described in subsection 2.2.2 were employed, which, considering the good agreement achieved with experimental results in section 5, appear sufficiently accurate for the present study.
The following properties, determined from the mesh technical drawing in Figure 6C, technical data sheets (Geobrugg, 2020)  Although the material behavior is in principle non-linear, the large deformations are predominantly rigid and accompanied by small strains, justifying the assumption of linear elastic material behavior.
By comparing Figures 8A,B a perfect agreement between the reference computer model and the real structure cannot be achieved as some posts in Figure 8B are already damaged. Additionally, the position and alignment of the supporting structure is unlikely to exactly match construction plans.

Impacting Object
In contrast to preceding works, such as Volkwein (2004), Volkwein (2005), Mentani et al. (2018) this work follows Sautter et al. (2021) to flexibly model arbitrary objects with discrete spherical element clusters. The advantage compared to standard finite element discretized objects is the simplified contact detection between arbitrary boundary objects and single spheres contained in the cluster as described by Santasusana et al. (2016). Figure 9 presents the DEM cluster models used within this work, and Figure 6B shows the real impacting object. While the clusters in Figures 9A,B represent the test object T092, Figures 9C,D show the cluster approximation of the test object T089. A study on the proper refinement level of such clusters can be found in Sautter et al. (2021). Special algorithms are required to create the refined cluster file, with this work employing the algorithms described in Bradshaw and O'Sullivan (2002), Bradshaw and O'Sullivan (2004) available in an online toolkit [http://isg.cs.tcd.ie/spheretree/].
Within this work several DEM particle properties are varied and their influence on the simulation result is studied, namely the friction, coefficient of restitution, and Young's Modulus. Contrasting this, the following physical properties are constant throughout all simulation runs: - Comparing Figures 6B, 9A it can be observed that the testing objects have already suffered damage from previous experiments. When comparing the simulation results with the ones obtained by the experiment in section 5 this unavoidable difference should be considered.
An explicit central-difference scheme as described by Matuttis and Chen (2014) was used to conduct the numerical time integration of the translational velocity of the cluster with a time step of dt 5 · 10 − 5 s. As presented in Irazábal et al. (2019), a more sophisticated time integration approach was used to integrate the rotational velocity using quaternions (Hamilton. 1844), which is especially critical due to the non-uniform geometry of the impacting clusters as shown in Figure 9.

VALIDATION
The following section presents the simulated Geobrugg 2017 experiment results to ascertain the practical applicability and accuracy of the aforementioned numerical modeling approaches and technologies. For clarity, key assumptions and uncertainties discussed in the preceding sections are summarized below: -The exact impact location was estimated from video records. -The model of the protective structure is taken from construction plans. However, it can be seen from photos and video recordings that the actual structure is partially crooked and already shows damage. -The reference angular orientation α DE of the impacting object describes the rotation around the main axis of the impacting object at the time of impact. It cannot be determined and is therefore included in the following investigation as an unknown variable. This also means that the angular velocities cannot be absolutely compared. The orientation of α DE is visualized in Figure 9A. -Three further DEM specific parameters cannot be taken unambiguously from the test and have to be varied to study their influence. Young's Modulus E DE , coefficient of restitution ϵ DE and friction μ DE are partly problemdependent and their influence is not well researched with regard to Attenuator barriers. -As visualized in Figure 6B the impacting test object already shows damage, which will be neglected in the modeling of the DEM sphere cluster.
The rock trajectory and the general trends of angular and translational velocities are used for comparison with the experiment.

Trajectories
The most reliable and functionally-important comparison is the simulated versus experimental rock trajectory. Attenuator barriers are used to protect exposed areas such as motorways from falling rocks, therefore the path of the object is of particular interest. In the following investigation the unknown model parameters were varied to check their influence on the results of the test run T092. To validate the results the test T089 was simulated subsequently with the best suitable parameters obtained from the investigation of T092. Figures 10, 11 illustrate the influence of the DEM parameters for T092  While the deformation of the impacting object is not of interest for this study the so-called Young's Modulus E DE represents only an algorithmic parameter in the calculation of the DEM contact forces (see subsection 2.1). The range of E DE in this work does not represent the physical properties of concrete but proved to result in the best fitting results. This could also be observed when the obtained parameters were used to simulate the test run T089 (see Figure 12B). Figure 11A demonstrates that E DE does not heavily influence this kind of simulation and further studies showed that choosing a E DE to be near the physical value of the Young's Modulus of concrete does not influence the simulation result but will lead to the necessity of much smaller times steps, contradicting the purpose of this work to offer a fast, efficient, simplified model and simulation of the Attenuator barriers.
An initial observation is that all variations give very good agreement with the test results for T092 and accurately predict the final position of the rock with a maximum error of +10%. While varying the reference angular orientation seems to have the least influence (see Figure 11B) the proper choice of a suitable friction value strongly influences the simulation as shown in Figure 10B. This sensitivity is expected as this experiment heavily   Figure 11A as well as the variation of the coefficient of restitution in Figure 10A show little influence, likely since the Attenuator barrier primarily retards via kinetic energy instead of strain energy.
In the following, only the three best parameter variants for T092 will be discussed, as the quality of the trajectories allows direct conclusions to be drawn about the correctness of the respective velocities. The three optimal parameter combinations are presented in the following list and their trajectories are visualized in Figure 12A together with the path of a free falling object to demonstrate the efficacy of the Attenuator barrier. The data obtained from the simulations is not only in strong accordance with the experiment results but also clearly demonstrates the effectiveness of the herein presented Attenuator barrier. A rock at height of 12m with an initial vertical velocity in gravity direction of 12.08m/s (as it is the case for T092) would need, 12.00m 12.08 m s · t + 1 2 · 9.81 m s 2 0.76s to reach the ground. In combination with an initial horizontal velocity of 10.01m/s the stone will have traveled a horizontal distance of 10.01 m s · 0.76s ≈ 7.6m. This results in a potential risk area that is approximately three times as large as for the protective installation (compare Figure 12A).
The results, which are shown in Figure 12B relate to the simulation of T089 with the parameter set C of T092. It can be clearly seen that the experimental results are properly reproduced without any further need of parameter investigation. This ensures that this parameter set can be used to recalculate further scenarios and creates confidence in the accuracy and applicability of the study presented here.

Angular Velocity
The experimental observations presented in Figure 13A illustrate the impacting object angularly decelerating and subsequently accelerating in the opposing angular direction. Unfortunately, the experimental data do not provide any information about predominate rotation axis, the individual rotation components of each axis, or the rotational orientation of the specimen at the time of impact. Nevertheless, it is useful to compare the general rotation behavior with the simulation results, with Figure 13B illustrating the simulation exhibits the same rotational trend as the experiment in Figure 13A. Additionally, Figure 14A visualizes the angular velocity about one of the main axes of the impacting object, which clearly illustrates the time period in which the aforementioned deceleration and subsequent acceleration take place. The same behavior can be observed in Figure 14B, which depicts the general trend of the angular velocity obtained from the experiment.

Translational Velocity
The translational velocity of the three optimal parameter sets (as per subsection 5.1) and the experiment are illustrated in Figure 15A. It can be seen that the deceleration of the impacting rock in the horizontal direction (perpendicular to the protection net) is accurately modeled by the simulation. This is of particular importance for the design of protection structures, since the management of horizontal momentum is the primary arresting mechanism. Although experimental results for the vertical (gravity) direction are substantially scattered, the simulation demonstrates a good agreement with the general trend of the experimental data.
Just as in subsection 5.1, the velocity components of the simulation and the experiment of T089 are also compared in the following. Figure 15B shows that these data could also be reproduced with very good agreement and thus allows suitable statements to be made about future rockfall events in Attenuator barriers.

CONCLUSION
This work has applied an advanced modeling and simulation approach to characterize the behavior of Attenuator barriers. Two independent numerical methods were coupled within an open-source multi-physics code to analyze the interaction of the impacting object and the protective structure. The impacting objects were simulated using DEM whereby arbitrarily shaped rocks were flexibly modeled as clusters of single spherical particles. This effective procedure allows the simultaneous use of simple contact algorithms and the possibility to consider the influence of arbitrary rock shapes. The structure itself was simulated using FEM employing a variety of advanced element formulations appropriate for each structural component. A homogeneous anisotropic membrane element was used to simplify the complicated protection mesh. As in Sautter et al. (2021), this simplification was found to be applicable when the global deformation behavior has to be investigated. A sophisticated yet elegant sliding cable formulation was applied to model the sliding nodes on the support ropes, allowing the boundary kinematics of the net to be fully resolved. Due to the robust implementation in the open-source code KRATOS [see Dadvand et al. (2010), Dadvand et al. (2013); Ferrándiz et al. (2020)], this work allows engineers and researchers alike to perform independent analyses of any rockfall event and the evaluation of suitable protection structures.
Another objective of this work was evaluating the performance of Attenuator barriers. It was demonstrated that this type of flexible protective structure offers dramatically increased protection against rockfalls by retarding impacts and subsequently guiding them to designated collection zones. In subsection 5.1 the simulation and experimental results agree that the Attenuator barrier considered reduced the danger zone by approximately two thirds compared to a free falling object. The fact that the Young's modulus and the coefficient of restitution of the DEM cluster have little influence on the trajectory of the impacting object was expected, since little energy is absorbed by strains and internal forces and the rock is instead deflected. On the other hand, the correct choice of interface friction was found to be important, as the experimental data is strongly influenced by the rotation of the impacting object. To support the conclusions of this study, another test was simulated without varying the unknown parameters. The set of best parameters from the previous tests was applied directly. The results of this work clearly show that with this procedure the experimental results can be calculated with a very good agreement. Thus, future questions regarding Attenuator barriers can be answered confidently, quickly and efficiently. The CPU system settings for this study is an Intel(R) Xeon(R) CPU E5-2,623 v4 at 2.60GHz, while it takes approximately 700 s (11.667 min) to run one simulation.
Live experiments are time-consuming and very costly, especially if every new change in the design has to be investigated experimentally. The numerical analysis of these tests can support the design process and thus make the development of new protective structures more cost-effective. Considering the global change of temperature and general environmental influences, the development and realization of simulation methods in open-source multi-physics codes (such as KRATOS) is of great importance and cannot be overemphasized. A community of engineers and programmers working together on such computer programs that are freely accessible to everyone is of paramount importance.

DATA AVAILABILITY STATEMENT
The data analyzed in this study is subject to the following licenses/restrictions: The comparative data in this work are the results of several experiments conducted by Geobrugg. The experiments were carried out in 2017 in Nicolum Quarry, British Columbia, Canada, with Wyllie and Norrish Rock Engineers and Geobrugg North America Ltd. The data are the property of Geobrugg and are partly published in the course of this work, the remaining data are not publicly available.

AUTHOR CONTRIBUTIONS
KS, PB, and RW have jointly developed the basic theory of the paper. KS has worked on the modeling, development and simulation of the multi-physics simulation discussed here. HH wrote the detailed explanation of the experiments and played a leading role in the ongoing discussions. CW, RW and K-UB provided the initial ideas and advised while developing the methods. KS, HH, CW and PW were involved in preparing and finalizing the paper.

FUNDING
This work was supported by the Technical University of Munich (TUM).