Partitioned Strong Coupling of Discrete Elements with Large Deformation Structural Finite Elements to Model Impact on Highly Flexible Tension Structures

(is article presents a staggered approach to couple the interaction of very flexible tension structures with large deformations, described with the finite element method (FEM), and objects undergoing large, complex, and arbitrary motions discretized with particle methods, in this case the discrete element method (DEM). (e quantitative solution quality and convergence rate of this partitioned approach is highly time step dependent. (us, the strong coupling approach is presented here, where the convergence is achieved in an iterative manner within each time step.(is approach helps increase the time step size significantly, decreases the overall computational costs, and improves the numerical stability. Moreover, the proposed algorithm enables the application of two independent, standalone codes for simulating DEM and structural FEM as blackbox solvers. Systematic evaluations of the newly proposed iterative coupling scheme with respect to accuracy, robustness, and efficiency as well as cross comparisons between strong and weak FEM-DEM coupling approaches are performed. Additionally, the approach is validated against the rest position of an impacting object, and further examples with objects impacting highly flexible protection structures are presented. Here, the protection nets are described with nonlinear structural finite elements and the impacting objects as DEM elements. To allow the interested reader to independently reproduce the results, detailed code and algorithm descriptions are included in the appendix.


Introduction
e simulation of lightweight structures impacted by heavy objects can be a challenging problem if numerical methods are used, as both sides have highly different masses. Some of the applications can be as follows: rockfall in protection net structures [1][2][3], racecourse protection structures to secure both spectators and drivers at the same time, suicide protection nets on bridges or towers, and many more. e examples in this work are mainly directed towards natural hazards, such as rockfall events which often cause destruction, especially in mountainous or populated areas.
As it is hard to prevent those events, protection structures are built along settlements and roads. In particular, lightweight flexible structures come into operation, as they are able to undergo large deformations and thus are capable of absorbing large amounts of energy with a smooth and comparably slow load transmission, reducing peak loads and maximum stresses. In principle, flexible structures are built to exploit the possibility to reduce stress peaks by their ability to tolerate large deformations. e interaction between the impacting object and very flexible structures calls for an advanced computational approach, as real scale tests are expensive, complex, and not suitable for a quick assessment of such structures. e abovementioned problems can be divided into two separate problems: the highly deforming fixed protection structures and the freely moving impact objects. Based upon both solutions, an interaction and equilibrium between both need to be found. e net structures can typically show high deflections. However, the topology usually does not change and thus the net can be described with finite element formulations having a meshed discretization. Furthermore, multiple kinematic formulations can be applied within the analyses. Approaches range, e.g., from shell or membrane structures [2,3] up to cable or beam models [4] or special formulations, such as ring elements [1]. In the following, the decision was made to use cable net models, discretized with cable element formulations (follow Section 2.2), as those provide accurate results for the given problems with comparatively small computational costs. Additionally, this approach allows contact at the correct position and thus allows small objects to pass through the structure.
Furthermore, the impact objects can also be discretized with different methodologies.
e current state-of-the-art rockfall impact simulations is to use rigid bodies to simulate impacting rocks on highly flexible structures [1] where damage and deformation of the impacting objects are neglected. e approaches include finite element methods (FEMs [1]), discrete element methods-(DEM-) described structures, or more flexible structures with material point methods (MPMs [5]), and others. Using the FEM makes it possible to accurately describe the continuum of impacting objects. However, the approach can be very costly, as the contact detection can become very complex. DEM provides much optimized contact detection and thus an efficient multiphysics simulation. e drawback is that the continuous expression is complex and dependent upon difficult parameter calibration. Furthermore, the use of MPM might show similar properties as the FEM with fewer contact difficulties. is method can be quite complex in terms of the setup and smears the contact due to its numerical properties. In the following, DEM is chosen to discretize the impact objects used herein, as an efficient algorithm is needed, and no further attention is focused on its continuum description.
Coupling the DEM and the FEM is a common way to simulate various multiphysics problems [6,7]. In particular, in problems where the interaction between granular materials or rigid objects with large motions and continuous systems is of interest, this combination of discretization methods is frequently used. Various applications, such as the thermomechanical behaviour [8] of contact between frictional bodies [9,10], assessment of strains in the simulation of shot peening [11,12], races and balls in ball bearings [13][14][15], general tribological systems [16,17] such as the simulation of rail tracks [18], and more advanced investigations, including fracture due to blast loads [19], are studied using DEM-FEM coupling. e DEM is also used to simulate production processes such as rotating machinery [10] and particulate flows [20]. e coupling between DEM and FEM is done in a partitioned manner, which allows the combination of the respective best-suited solution strategies for each subproblem and the transfer of information in-between.
Accordingly, the user is not restricted to a code which includes both participants but instead can couple any existing DEM and FEM software by creating a suitable interface. is publication concentrates on the discussion of spatial mapping with cable structures. e coupling ultimately also allows the use of blackbox solvers for each participant (e.g., symplectic Euler [21] for the DEM and generalized alpha method [22] for the FEM). In order to advance the state of the art, a strong coupling algorithm is developed for the DEM-FEM impact problems.
To investigate its potential and superior performance in the underlying application case, a weak coupling algorithm is presented and used for comparison. It is known from the field of fluid-structure interaction [23,24] that a strong coupling algorithm typically allows larger time steps compared to a weak coupling approach. e aim of this work is to investigate the properties of a strong coupling algorithm for DEM-FEM coupling and assess its usability in this setup.
From a formal point of view, the structure of the paper is as follows: (i) Section 2 describes the FEM notation including an introduction to the applied cable formulation (ii) Section 3 gives an overview of DEM (iii) Section 4 introduces the equilibrium required for the coupling and depicts the spatial mapping (iv) Section 5 explains the staggered weak coupling approach (v) Section 6 depicts the strong coupling approach, which adds additional complexity (vi) Section 7 demonstrates the advantages of the proposed coupling algorithm and investigates the influence of a variety of different input parameters (vii) Sections 7.1 and 7.2 show and investigate the novel coupling approaches (viii) Sections 7.3 and 7.4 present sensitivity analyses of the important parameters (ix) Sections 7.5 and 7.7 demonstrate large-scale problems (x) Section 8 gives a conclusion and outlook on future research (xi) Appendix A provides a more detailed discussion of the calculation of contact forces for the DEM (xii) Appendix B provides references to the software used in this study and provides scripts for reproduction of results

The Finite Element Method
e finite element method (FEM) [25] is used to numerically solve partial differential equations. For this purpose, a domain is discretized into finite elements for which an approximated solution is known.
As described in [22,[26][27][28], kinematic relations describe the possible movement of such elements. e current configuration x can be described with the help of a time-2 Advances in Civil Engineering dependent map (noted as t) from the reference configuration X, which can be obtained from the displacement vector u (additionally see Figure 1): the bold letters indicate tensors and vectors.

Virtual
Work. e momentum equation can be described by applying the divergence theorem [22,27,29,30]: with σ being the first Cauchy stress, the body forces b, the acceleration € u, and the density ρ. Multiplying the momentum equation (2) with the virtual displacement field δu and integrating over the physical domain Ω leads to the weak form, respectively: e three virtual work contributions: the internal work (equations (4), respectively (5)), the external work (equation (6)), and the kinetic work (equation (7)), are obtained by the partial integration: with the Almansi strain e, the normal direction of Ω, n, and the Neumann boundary Γ σ . From equation (4), one can obtain equation (5) (for detailed information, see [22,27]), which describes the internal virtual work with the use of the Piola-Kirchhoff 2 (PK2) stress tensor S and the Green--Lagrange strain tensor E within the initial domain Ω 0 .

Cable Element.
With regard to the cable element, a cable or a "truss element formulation" having two nodes (n k and n l ) considering geometrical nonlinear behaviour is used. e kinematics is outlined in Figure 1. Each node requires three (or for 2D, two) displacement degrees of freedom (dofs), but no further dofs are required. As the bases for the structural analysis and to introduce relevant notations, the internal forces and the tangent stiffness matrix are briefly noted.
Due to the one-dimensional nature of the element, only one parametric coordinate ξ [−1 ≤ ξ ≤ + 1] is used to describe the geometry [25] with linear shape functions.
By virtue of equation (5), the internal forces F int read for each degree of freedom r as follows: for the constant crosssection A, the reference length L, and a given prestress S pre . e tangent stiffness matrix K can be expressed as e constitutive law used in this work is the St. Venant-Kirchhoff hyperelastic material model [29,30] described by the strain-energy functional Ψ sv . Its first and second derivatives w.r.t. E describe the PK2-stress S and the material tangent modulus D, respectively [30]: e prestress S pre is a stress, which is initially applied and added to the internal forces of the structure. It is a secondorder tensor; however, for the truss kinematics, only a scalar is required. If a stress state, in equilibrium to the initial prestress, is required within a complex structure, formfinding procedures need to be applied beforehand. In contrast to standard structural analysis, form-finding solves an inverse engineering problem. is inverse problem is to find a geometry which corresponds to an equilibrium configuration, whose stress state is identical to the given prestress. e force density method is one approach used to realize this [24,31]. e described kinematic relation is capable to model either a cable or a truss formulation. To model a realistic behaviour of cable structures, it needs to be ensured that no compression stresses are carried by the structure. With an additional check if a cable is under compression, stiffness contributions and internal forces are avoided. In the following, using the applied naming, trusses would allow compression stresses, whereas cables do not. e element mass matrix M for a dynamic simulation results from the virtual kinetic work as depicted in equation (7). With N(ξ) being the shape functions and the reference density ρ 0 , M is defined as follows: e damping matrix C is expressed via the Rayleigh damping as a linear combination of the mass matrix M and the stiffness matrix K (for more information, see [32]): with τ and κ as the scaling factors.
Advances in Civil Engineering 3

The Discrete Element Method
is section describes the fundamentals of the discrete element method (DEM), with a focus towards the DEM-FEM coupling (following Sections 4-7). e DEM used within this scope models the dynamics of spherical particles, considering external forces such as gravity and contact with other particles or contact with differently described objects. DEM can be used for many different particle shapes such as rectangles, cones, spheres, and more. However, only spherical particles are considered in this publication. To describe more complex shapes, a set of spheres is connected within clusters. A cluster consists of multiple particles which are used for contact detection and force evaluation. e mass and centre of gravity are described within the cluster shape and independent of the masses of the particles (for further information, see [33]). e basic steps in a DEM simulation are as follows (also [21]): (i) Contact detection (ii) Evaluation of contact forces (iii) Integration of motion ese steps will be briefly described in the following.

Contact Detection.
Within the scope of this publication, certain contact scenarios are of interest. e centre of spheres is described with C i , and their respective radius is R i . e distance to contacted edges is deciphered with d i (Figure 2). Two spheres C i and C j : Sphere C i and vertex n k : Sphere C i and d i to edge:

Evaluation of Contact
Forces. e contact forces can be evaluated using different contact laws and rheological models. For a detailed description of contact models, see [6,21,34,35]. A more detailed description of the evaluation of contact forces is included in Appendix A. e coefficient of restitution (COR), ε ⊥ (see equation (A.7)), is an essential part in the contact and is further investigated in Section 7.4.

Integration of Motion.
e integration of motion is described by Newton's second law of motion. e translational acceleration € u is described via the force F and the mass of the particle m. e angular acceleration _ ω is expressed by the moments T and the tensor of rotational inertia I [21]: e forces and torques at each particle i are described as follows [6,19,21]: e resulting forces and torques depend on the following components: F ext i , T ext i : external loads (e.g., gravity) F ij : contact interactions between neighbouring particles or boundaries as the contact F contact (see Appendix A for description of the components F n and F t )  : external damping boundaries n ij , t ij : normal and tangential vectors at the respective contact point r ij c : vector connecting the particle centre and the respective contact point

Staggered Coupling of DEM and FEM
To couple two standalone physically independent interacting numerical problems such as the simulation of particles and structures, a suitable interface condition to achieve equilibrium is needed to be found.

Structure-Particle Equilibrium.
To put the two independent physics, the DEM and the FEM, into an equilibrium state, the following force equilibrium at the structure needs to be achieved: where u describes the displacements, _ u represents the velocities, Ω S represents the structural domain, while Γ of Ω S includes all nodes which are in interaction with the DEM particles, P.
With respect to equation (18), the contact forces F contact of the particles, which are dependent on its displacements and the velocities, need to be in equilibrium with the internal forces F Ω S,Γ int of the structure, which are dependent on its displacements, velocities, and accelerations (equation (8) and additional including inertia effects). e equilibrium of both DEM and FEM simulation is graphically depicted in Figure 3. e basic idea of the proposed partitioned coupling simulation is the interchange of primary (such as the displacement) and secondary (e.g., forces) interface variables which are obtained as the solution of the respective components of the simulation.

Spatial
Mapping. In the following, the DEM problem is solved independently from the structural problem. To do this, the displacements and velocities of the structure at the given time step are transferred to DEM, and this structure is further seen as the DEM wall, described by the domain Ω D . e wall is used to calculate contact forces with the DEM particles P, which are depending on its displacements and velocities (see Appendix A for contact laws and force calculations). After solving the DEM problem, the resulting contact forces are transferred to the structural analysis problem. With the contact forces, seen as external forces, the dynamic structural problem is solved, resulting in new displacements and velocities on the domain Ω S . is procedure is outlined in Figure 4.
Following this, the contact forces F contact are now dependent on the displacements and velocities of Ω D and not directly on Ω S and are defined as the external forces coming from the DEM: e equilibrium within the structural mechanics problem is given as follows: After solving both domains, the two interface conditions, for the displacements and the velocities between both fields, are not fulfilled anymore: Resulting from this, the contact forces F contact computed within Ω D and the contact forces which would be computed within Ω S are not the same anymore, and thus, the equilibrium expression is not fulfilled: For small time steps, resulting into smaller contact forces, the tracking of the interface equilibrium can be negligible. However, for ill-conditioned systems and large time steps, the resulting difference will lead to inaccuracies

Advances in Civil Engineering
and makes the solution unstable. To solve this problem, a possible approach is presented in Section 6.

Influence of Coefficient of Restitution (COR).
Large contact forces will result in difficult fulfilment of interface conditions (equations (21)-(23)). Section 6 proposes a remedy for that problem. One major factor influencing the magnitude of the contact forces is the DEM particle property COR. is value must be defined by the user and heavily influences the stability of the coupled simulation (see example in Section 7.4). e coefficient represents the ratio of initial speed and final speed after impact [21] (equation (A.6)) and is further discussed in Appendix A. Since this coefficient is determined manually for each simulation, it is important to be careful when doing the calibration.

Mesh Dependency for Cable Structures.
For the specific application of highly flexible cable structures in this study, such as rockfall protection nets or any other kind of cablelike structures, the DEM wall condition Ω D discretization and the FEM Ω S discretization on Γ must exactly coincide (conforming meshes). e respective meshes represent a physical mesh which must be correctly described in order to model the exact contact positions. To demonstrate this behaviour, Figure 5(a) visualizes the DEM part of the simulation. A cable net is modelled and impacted by two spheres. A large sphere finds contact and deforms the boundary, while a smaller sphere penetrates an opening. In addition, Figure 5(b) represents the respective FEM structure which is used to calculate the adequate structural response.
If surface elements such as shells or membranes, which do not possess physically the predetermined discrete contact positions, are used within a coupled simulation, arbitrary meshes can be used. In that case, a mapper [36] will be responsible for the correct data transfer.

Staggered Weak Coupling
e fundamental idea of the weak coupling (sometimes also called explicit coupling [24]) follows a single exchange of coupling data in each time step. e communication pattern is depicted in Figure 6. e important steps at each time, including this communication pattern, can be summarized as follows: (1) Solve DEM (results: u P , _ u P , and F contact ) (2) Map F contact from DEM to Structure (3) Solve Structure (results: u Ω S,Γ , _ u Ω S,Γ , and € u Ω S,Γ ) (4) Map displacements and velocities from Structure to DEM (5) Advance in time (not explicitly shown) e interface variables are accordingly updated (see Steps (2) and (4)): Velocity: Contact force: is algorithm is comparatively easy to implement and typically does not require any deep interaction. Standard DEM and FEM simulation environments provide the exchange data as an output. erefore, different software can also be efficiently applied here. Furthermore, it was shown that the algorithm can be applied if the time steps do not become too large (see examples in Sections 7.1 and 7.2). However, the behaviour of this procedure can become unstable as soon as the differences in stiffness, mass, and velocity between the two physics become very high. e procedure is then very prone to the time step size used. Decreasing the time step size will lead to inefficient and numerically costly simulations.
To gain a deeper understanding of the underlying procedure, this approach is further detailed in Algorithm 1. In order to facilitate the reproduction of the results, the Python script used, including all comments, is provided in Appendix B.
In this procedure, two additional features will be discussed. ey are independent of the coupling approach used but improve the performance significantly. ey are added in Algorithm 1 and highlighted as follows: (i) particle_near_wall (line 3-6): if the respective particles are in the vicinity of the structural model to adjust the time step is checked. A particle moving freely in space can be simulated with a time step larger than it would be required for the simulation of the DEM-FEM interaction. (ii) forces ≠0.0 (line [12][13][14][15]: this is an additional check used only to solve the structure when contact forces are present. is is only valid if the self-weight of the structure or any other loads (except impact loads) are neglected. Otherwise, a preliminary simulation or a form finding [24,31] of the structure is needed.

Staggered Strong Coupling
As known from other coupled multiphysics problems, such as fluid-structure interaction (FSI) [24,37], the direct explicit transfer of the interface data (forces, velocities, and displacements) can lead to divergence problems in the staggered simulation. is problem is caused by large contact forces due to differences in velocities, acceleration, and highly different masses on both sides. In contrast to the weak coupling approach, the strong coupling (in the literature also called implicit coupling [24] or a conventional serial staggered approach within the context of loose coupling [38,39]) adds an additional iteration loop in each time step, which solves for the equilibrium between both numerical physics. is requires a Gauss-Seidel loop between DEM and FEM, which might need to be solved multiple Advances in Civil Engineering times within one time step [23,24,38,39]. is strategy enforces the coupling conditions (equations (4)- (6)) to be fulfilled. Convergence is considered to be achieved, as soon as the interface residual is below a user-defined tolerance ε. e residual formulation is defined by using equation (27) e steps of this approach are shown in Figure 7 and summarized in the following, using the respective numbering in the abovementioned Figure 7: (1) Solve DEM (results: u P , _ u P , and F contact ) (2) Map F contact from DEM to Structure (3) Solve Structure (results: u Ω S,Γ , _ u Ω S,Γ , and € u Ω S,Γ ) (4) Map displacements and velocities from Structure to DEM (5) Calculate interface residual ϵ (equation (27)) (6) Repeat Steps 1-5 until the interface residual reaches a given tolerance (7) Advance in time e weak coupling algorithm, described in the preceding Section 5 expresses single iteration in the strong coupling scheme (Steps (1)-(4)). e additional interface loop (Step (6), being controlled by the breaking criteria in Step (5)) which adds complexity to the solution procedure and significantly increases the computation costs as the system now needs to be solved multiple times within one time step. However, it allows more accurate results and higher simulation stability. It can be noted that the number of solving iterations is typically still lower than if the time step would be reduced to a value where the weak coupling approach would still be applicable. is is especially due to the property that  (1) Initialize Solve structure (FEM) (14) Map velocity and displacement on Γ from Ω S to Ω D (15) Update position of Ω D to Ω S ⟶ equation (1) (16) Finalize ALGORITHM 1: Weak coupling. 8 Advances in Civil Engineering many coupling iterations are typically not required throughout the simulation, but only at specific time steps. e comparison of the two procedures, including a view on the performance, is outlined in Section 7. e residual criteria within the strong coupling loop are defined by where ϵ is the user-defined breaking tolerance. It is checked after each iteration k by scaling the norm of the residuum k r with the square root of the number of degrees of freedom n eq at the interface Γ [40]. It is important to note that the interface tolerance should be larger than the convergence tolerances of the respective individual solvers within the coupled system; otherwise the convergence criteria cannot be reached. e residuum can either be obtained by the displacements, the velocities, or the contact forces. By subtracting the current solutions on the boundary Γ from the previous solutions of Step k − 1, the residuum of each variable can be noted as follows: Velocity residuum: Contact force residuum: Furthermore, large time steps typically lead to large differences in the interface velocities and displacements, and thus the result can be nonphysical large contact forces. If those forces are too high, small time steps still can lead to unstable simulations, even with the use of the proposed strong coupling algorithm. As a remedy, the transferred data can be gradually applied, which is also called relaxation. e outcome is that this permits a faster interface convergence. e so-called convergence acceleration [24] can be achieved by numerous methods and is discussed in the following.
Two different strategies can be chosen for the relaxation: either the relaxation of the displacements and velocities or the relaxation of the contact forces. e relaxation is done w.r.t. the residual (equations (28)-(30)), respectively: Relaxed velocities: Relaxed contact forces: Each variable is subsequently updated from the previous solution (Step k − 1) using the respective interface residuum scaled by relaxation factor α.
ere are different approaches to obtain the scaling factor α [39]. e relaxation factor can be set to a userdefined constant value, which is very simple and helps to improve the quality of the simulation. Another approach is to use the Aitken method [41]. e Aitken method optimizes α in every iteration w.r.t. the current residuum k r and the previous residuum k−1 r: respectively, α u � α(r u ), α _ u � α(r _ u ), and α F � α(r F ). e influence of the relaxation factor α is studied in the example in Section 7.3.
In this study, either the displacement and the velocity field or the contact forces are independently relaxed and subsequently mapped. However, in the case of displacements and velocities, both residua have to be achieved to ensure that both solution fields still coincide on both sides. us, the resulting residuum for both relaxing procedures is given as follows: Displacement and velocity residua: Contact force residuum: e interface variables are updated accordingly (see Steps 2 and 4 in Figure 7). e following variables are exchanged within the interface:

Without relaxation
Displacements: Velocities: Contact forces: With relaxation Displacements: Velocities: Contact forces: In summary, both solution strategies are described within Algorithms 2 and 3 in pseudocode. ey are both further elaborated on in Appendix B.3.

Systematic Assessment of the DEM-FEM Coupling
is section presents some examples which systematically analyse the difference between the herein introduced coupling approaches and their application within the simulation of relevant industrial applications.
e examples show problems of impacting objects on highly flexible lightweight cable structures, such as protection nets. ese interaction problems typically have numerical stability issues within the simulations, as the net structures have a low mass, whereas the rocks are typically heavy. is instability leads to the problem that especially when the first impact occurs, the forces might become very large. us, due to the different masses, this may lead to convergence problems, especially if the chosen time step is large, which can lead to inaccuracies in the simulation.
In the first academic problem 7.1, a cable structure is modelled to evaluate the influence of different time step values. Section 7.2 subsequently uses a cable structure with a large prestress while also showing the influence of the COR in order to analyse the influence of larger contact forces on the required time step. Section 7.2 investigates the difference between relaxing forces (Algorithm 3) and relaxing displacements and velocities (Algorithm 2). e proper choice of a relaxation factor is further discussed in the example in 7.3. e influence of the COR, which scales the contact forces, is then analysed in Section 7.4. Finally, a practical application of a rockfall into a cable net, using the herein explained approaches, is presented in Section 7.5.

Impact on a Compliant Cable: Large Deformations.
In this example, a single DEM particle with perfect spherical dimensions impacts on a prestressed cable, which is discretized with three finite elements. Here, the contact point on the structure is known, and thus it can be focused on the performance of the coupling algorithms. e setup of this academic example can be found in Figure 8(a), with E as Young's Modulus and Poisson's ratio ]. It demonstrates the necessity of a strong coupling algorithm; since for larger time steps, the phenomena of artificial contact loss, due to large initial contact forces, occur.
Within empirical tests, the time step Δt � 10 − 3 s is found to be the highest possible time step for which the weak coupling algorithm can resolve to an appropriate solution.
(1) Initialize (2) While time < end_time do (3) While interface _ res ≥ tolerance _ interface do (4) Search nearest neighbours and find contact ⟶ equations (13)-(15) (5) Calculate contact forces ⟶ equation (A.5) (6) Time integration of DEM part ⟶ equation (16)  (7) Map forces on Γ from Ω D to Ω S (8) Solve structure (FEM) (9) Map velocity and displacement on Γ from Ω S to Ω D (10) Calculate residuals for velocity and displacement ⟶ equations (28) and (29)  (11) Relax velocity and displacement ⟶ equations (31) and (32)  (12) Update position of Ω D ⟶ equation (1)  (13) interface_res � max(displacement_residual, velocity_residual) Time integration of DEM part ⟶ equation (16)  (7) Map forces on Γ from Ω D to Ω S (8) Calculate residuals for forces ⟶ equation (30)  (9) Relax forces ⟶ equation (33)  (10) Solve structure (FEM) (11) Map velocity and displacement on Γ from Ω s to Ω D (12) Update position of Ω D ⟶ equation (1) (13) interface_res � force_residual Here, the coefficient of restitution (COR) is set to be 1.0. Implicit time integration is used for the structure, as the chosen time step is too large for an appropriate solution with an explicit time integration scheme. Figures 8(b) and 8(c) show the behaviour of the cable after impact, for the time step size of Δt � 10 − 2 s. Using the weak coupling approach, a "jump" can be outlined as shown in Figure 8(b). Due to the large time step, a greater indentation and higher velocities occur. Consequently, the interacting force is too large, so that the sphere and cable do not continuously stay in contact during the entire time. is leads to a nonphysical behaviour of the coupled problem, as shown in Figure 9. By adding the additional interface loop to solve for the contact force, the convergence of the problem can be achieved for a larger time step of Δt � 10 − 2 s.
In the following, the time step of the first contact is discussed in detail. It can be seen (Figure 10(d)) that the contact force is relatively large in the first inner iteration (coming from the relatively large time step) and decreases within the interface iteration to a converged solution, due to the application of the Aitken relaxation, introduced in equation (34). is exemplarily demonstrates the advantages of the strong coupling scheme, presented in this article. e large discrepancy in the contact force would lead to an unstable coupled simulation when using a standard weak coupling algorithm. e same accounts for the deflections of the impacting sphere as shown in Figures 10(a)-10(c) presenting a visual description of the interface condition in equations (21) and (22) It can be seen that the positions of Ω D and Ω S do converge to a common value to fulfil the interface displacement/velocity equilibrium (equations (21) and (22)).
As an example, the converging contact force for each iteration in time step t � 1.2 s and the total number of inner iterations for each time step are shown in Figures 11(a) and 11(b). It shows that the force at the first iteration is almost 4 times higher than it is after the converged solution. Figure 11(b) shows that the number of inner iterations can vary greatly (between 1, if there is no contact, and 9 iterations) within each time step. However, it can be noted that the number of contact simulations is still lower than if the time step would be decreased to Δt � 10 − 3 s, which is the limit for which the weak coupling approach still converges.

Comparison to Position of Rest with Different Time Steps.
In this section, a setup similar to the previous example (Section 7.1 and Figure 8(a)) with an increased prestress (S pre � 10 6 N/m 2 ) and a homogeneous Young's Modulus (E � 10 9 N/m 2 ) in the cable structure is used with the following changes for the impacting sphere: R � 0.12 m and ρ � 3.5 · 10 4 kg/m 3 . e result of the transient analysis will be compared to the static solution, considering the particle as an external force. is static force is defined as follows: (43) e resulting static deflection of Point A (Figure 8(a)) is shown in Figure 12.
is comparison proves that the transient analysis approaches the static solution after a certain time.
Furthermore, the sensitivity of the time step within each coupling algorithm is also studied in this example. e results of all solutions are presented and compared in Figure 12. It shows that the weak coupling approach provides an accurate performance for a time step of Δt � 10 − 3 s, whereas the solution for Δt � 10 − 2 s is very unstable. It turns out that for large time steps, the result oscillates around the expected solutions. e measured solutions for time steps of Δt � 10 − 2 s and Δt � 3 · 10 − 2 s show that the strong coupling algorithm still allows for good convergence for rather large time steps. However, by increasing the time steps, the number of interface iterations subsequently increase, which is shown in Figures 13(a) and 13(b). Especially when the impact occurs, the large difference in interface velocities leads to an increased number of interface iterations (Figure 13(b)). Within the scope of coupled simulations, this difference is decreased by the proposed algorithm leading to a smaller number of iterations. e influence if either displacements and velocities or forces are relaxed is examined in the following. Both options are described in Algorithms 2 and 3, respectively.
Comparing Figures 13(a) and 13(b), it can be noted that relaxing the forces facilitates slightly faster convergence than relaxing displacements and velocities.
In this specific case at hand, clear and marked off points of load application (impact position) do exist. In different cases, for example in the following Section 7.5, where a variety of possible impact nodes exist, relaxing displacements and velocities are shown to be the better choice. In those cases, which appear more frequently, the impacting spheres can rapidly change the impacting position and thus lead to a slow converging force residual.

Influence of the Relaxation Factor.
In this example, the influence of the relaxation factor α in the case of relaxed displacements and velocities is investigated by comparing the Aitken relaxation (equation (34)) and a set of constant relaxation factors. For this purpose, the COR is set to 1.0 (which physically describes a perfectly elastic impact on a rigid obstacle) to neglect the influence of the wall velocity and thus concentrate solely on the relaxation factor. e same problem setup as in Section 7.2 is used.
As Figure 14 shows, a constant relaxation factor can be used as long as it is smaller than 1.0. α � 1.0 describes a nonrelaxed system and does not find a proper solution for 12 Advances in Civil Engineering Advances in Civil Engineering 13 this given example. Manually finding a suitable constant relaxation factor is cumbersome and is dependent on the system setup. In addition, it heavily influences the solving time, as Table 1 demonstrates. For a constant α and the Aitken α (equation (34)), the comparison is performed with respect to computation time. It can be noted that although constant relaxation factors provide good results, the optimized Aitken relaxation factor facilitates faster convergence to the residual.

Influence of the Coefficient of Restitution.
Another important entity within the multiphysics problem is the COR ε ⊥ (Section 3). As Figure 15(b) shows, the COR directly influences the contribution of the DEM rigid wall velocity to the contact force. Current state-of-the-art publications such as [42][43][44] express the importance of the COR value for impact simulations. For a case study, different COR values are used, while the time step is kept constant. As can be seen in Figure 16(a) (zoomed in Figure 16(b)), the interface coupling becomes unstable as soon as the COR reaches a small value. is instability can be overcome by using the strong coupling algorithm presented in Section 6 and is a result of the increased contact force in the system [34]. Additionally, Figure 16(a) describes another important feature: the choice of COR does not influence the final damped solution of the structure (see "static" in the graph in Figure 12) but only the maximum transient solution.   Figure 12: Comparison of weak coupling approach with Δt � 10 − 2 s and Δt � 10 − 3 s against the strong coupling approach with Δt � 10 − 2 s and Δt � 3 · 10 − 2 s for the vertical displacements of node A (Figure 8(a)). Furthermore, the static solution of the problem is provided to show which value the results should approach. additionally introduced interface iterations. As the simulation proceeds and the initial velocity difference is properly controlled, fewer iterations are needed to enforce the interface conditions.

Practical Application: Angled Protection Net.
One prominent practical application case of highly flexible structures can be found in mountainous regions. As an  With α � 1.0, the simulation becomes unstable shortly after the first contact is detected. Underrelaxation of α < 1.0 leads to stable simulations too. However, this is more computationally costly than using the optimized Aitken relaxation factor.  alternative to protection nets used to catch rocks, angled nets can also be spanned over roads to direct impacting objects to a safe spot, as shown in Figure 18(a).
To test the limits of the presented algorithms, in this study, the same system as shown in Figure 18(b) is modelled without prestressing the cable structure, leading to a very compliant structure (compare Table 2). Additionally, a small COR of ε ⊥ � 0.2 and a high impact velocity are chosen to introduce even more difficulties due to an increased contact force.
Using a time step of 2 · 10 − 4 s, the different behaviours after impact are presented in the incidental Figures 19(a) and  19(b).
Similar to the example from Section 7.1, the weakly coupled problem experiences too large contact forces and loses contact between the impacting object and the structure, whereas the strong coupling algorithm manages to keep the contact for the given time step (Figure 20(a)). e considerably large number of interface iterations (Figure 20(b)) is a result of the system setup. is example Advances in Civil Engineering tries to push the time step to a maximum and represents the largest possible time step, which cannot even be used for weak coupling anymore, describing a complex problem.

Special Modelling Possibilities.
Using two standalone solution techniques, such as the DEM and structural mechanics FEM, enables the user to benefit from the full range of capabilities and strengths of both participants, such as sliding nodes on cable elements (including friction) [1,4] ( Figures 23(a)-23(c)), custom ring elements [1,4,46], plasticity laws to model energy dissipation elements, choice of multiple time integration schemes, and clusters of particles (Figures 24(a)-24(d)) to model arbitrarily shaped objects such as rocks, which is an advantage over state-ofthe-art rockfall protection simulations as discussed in [1,47].
Using rigid bodies to model impacting objects has the disadvantage of neglecting damage and deformation of the object itself and adding additional complexity when handling arbitrarily shaped objects. e DEM offers the possibility to simulate breakup of impacting objects [21]; however, the simulation of the continuum with DEM particles can be very costly, and calibration can be elaborate. e possibility to model Ω D with line wall conditions suitable to the FEM mesh additionally allows for small rocks to penetrate the protection net (Figures 25(a)-25(e)). In summary, the combination of the DEM and the FEM allows the user to model any possible object impacting in a highly flexible structure modelled by any suitable structural finite element types. Any conceivable combination ( Figure 26) is possible as long as a suitable algorithm, as presented in this paper, is available.

Conclusions and Outlook
e numerical analysis of lightweight structures coupled with impacting heavy objects proves to be a complex problem and leads to instabilities within the simulation, especially due to the different masses of the participants. To overcome this problem, this publication presents several staggered coupling approaches and presents a sensitivity study with respect to certain crucial parameters. e procedure suggested herein uses FEM with cable element formulations for flexible lightweight structures (Section 2) and DEM for the interacting objects (Section 3). Furthermore, Section 4 shows the procedure for reaching the equilibrium between both physics. First, the procedure is explained with a single interface calculation within each time step (Section 5). In many examples, this approach proved to be unstable, specifically at initial contacts (indicating large velocity differences). Additionally, the simulation needs small time steps, which might be required only at certain steps. us, to overcome this problem, an additional iteration between the physics was explained in Section 6. is allows time steps to be increased significantly and improves the efficiency of the simulation (see examples in 7.1, 7.2, and 7.3). Additionally, the sensitivity of the quality of the simulation is tested by varying the relaxation factor (equation (34) in combination with example in 7.3) and the coefficient of restitution (COR) (see Appendix A in combination with example in 7.4). While the underlying algorithms are abstractly presented in the preceding sections, more detailed versions can be found in the following appendix to allow the interested reader to independently reproduce the results. e novel approaches make it possible to efficiently simulate the correct behaviour of complex existing structures. e example in 7.5 shows net structures interacted with rocks which are based on existing structures in the Austrian and Swiss Alps. e stability is heavily influenced by a restricting time step (Figures 19(a) and 20(a)) if the interface is not controlled by a suitable algorithm as presented in this study.
In addition, the use of two standalone applications in this study, the so-called "blackbox solvers" allows for a variety of advantageous features. As described in Section 7.6, any given terrain model can be integrated into the simulation process to efficiently capture environmental influences on the results (Figures 22(a)-22(e)). Furthermore, Section 7.7 demonstrates the advantages of an independent FEM application which is capable of modelling numerous structural details, such as energy dissipation elements or sliding nodes on cable elements. Accordingly, DEM can be used to model arbitrarily shaped impacting objects (Figures 24(a)-24(d)). is allows for independent work in the respective application without changing the coupling strategy, which especially proves beneficial in an open-source software environment such as KRATOS [48].
In future research, different FEM formulations [2,3] can be tested for the simulation of the protection nets. Furthermore, if rocks cannot be explicitly described, other particle approaches such as the material point method (MPM [5]) could be applied with the proposed coupling approach. By the way of example, conceivable application cases include the simulation of mud-flow barriers as well as avalanche barriers. Furthermore, the influence of the time integration scheme is additionally a significant factor which will require deeper investigations. Fluid-structure interaction.

Abbreviations
In order to obtain Δs, the tangent unit vector t ij must be derived by splitting the velocity at the contact point v and then, split which allows us to express Δs using the tangent unit vector, the element displacement u i , and the element rotation Θ i : For a Hertz-Mindlin spring-dashpot contact model (denominated as HM + D in [34]), as shown in Figures 27(a) and 27(b), the normal F n and tangential F t contact forces [21,34] for the case of two spheres colliding, F n � k n δ n + c n _ δ n , F n te � F n−1 te + k n t Δs n , if ΔF n ≥ 0, are calculated with the aid of the normal and tangential stiffness k n and k t , respectively, and the damping coefficients c n and c t [21], considering the maximum tangential force restricted to the Coulomb's friction limit [34] with the coefficient of friction μ. e tangential forces are consequently updated from the last step, indicated by the superscript n. e material parameters in Table 3 (Young's modulus E, particle mass m, shear modulus G, and Poisson's ratio ]) are typically obtained by calibration from experiments [49].
With respect to [21,35], the COR expresses the ratio between the velocity after _ δ after n and the velocity before _ δ before n impact. A maximum of ε ⊥ � 1.0 will model a perfectly elastic impact, whereas ε ⊥ � 0.0 models a perfectly plastic impact. A smaller COR increases the influence of the impact velocity in the contact force calculation [21,34] (equation (A.5)) and is thus critical for the coupled simulation in this study, in which a Hertz-viscous-Coulomb contact law is used [21]. For frictional cohesion-less contact as used in this study, the normal force must be constrained to always be ≥0.0 [21], since no traction in normal direction is allowed. To correctly capture this behaviour, Cummins et al. [34] describe how to modify the contact duration calculation by using the dashpot coefficient as described by equation (A.8). For further discussion on the contact forces, additional contact laws, and specific DEM-related topics such as contact duration, etc., see [21,34,35,50,51].

B. Code Scripts and Development Environment
To give the interested reader a better understanding and the possibility to reproduce the results, the algorithms are presented with the notation used. e open-source multiphysics software KRATOS [48] was used for this study. It can be downloaded [52]. An installation guideline is provided there, too. KRATOS is designed in C++ and includes a Python interface to facilitate the advanced development and simulation. Documentation for the Python scripts used in this study is provided in the following. To run the simulation, the structural mechanics application, the discrete element application, and the mapping application are required.

B.1. Problem Setup.
e script to define the problem setup is shown in the following. is initialization script is for both the weak and the strong coupling approach, which are described in the following appendices (Algorithm 4).

B.2. Weak Coupling Algorithm.
First, the Python script to run the weak coupling algorithm (Section 5) is provided.
is code sequence describes two possibilities to improve the efficiency of the simulation: one by increasing the time step if particles are far away from the interface and the other by solving the FEM part only if contact forces have been detected. ese two features are omitted in Appendix B.3 for simplicity purposes. However, they can be used to optimize computation time (Algorithm 5).

B.3. Strong Coupling Algorithm.
e two strong coupling approaches described in Section 6 are depicted in the following algorithms. First, the procedure to relax displacements and velocities is explained, followed by the procedure to relax the transferred forces (Algorithms 6 and 7).

Data Availability
No data were used to support the findings of the study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.

Authors' Contributions
All the authors prepared the manuscript. All the authors read and approved the final manuscript.