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

Authors: Klaus Bernd Sautter; Tobias Teschemacher; Miguel Ángel Celigueta; Philipp Bucher; Kai-Uwe Bletzinger; Roland Wüchner
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

This 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). The quantitative solution quality and convergence rate of this partitioned approach is highly time step dependent. Thus, the strong coupling approach is presented here, where the convergence is achieved in an iterative manner within each time step. This 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.


---

## Body

## 1. Introduction

The 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–3], racecourse protection structures to secure both spectators and drivers at the same time, suicide protection nets on bridges or towers, and many more.The 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. The 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.The 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. The 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. The 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. The 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. The 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. This 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–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. The DEM is also used to simulate production processes such as rotating machinery [10] and particulate flows [20].The 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. This publication concentrates on the discussion of spatial mapping with cable structures. The 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. The 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)
Sections7.1 and 7.2 show and investigate the novel coupling approaches(viii)
Sections7.3 and 7.4 present sensitivity analyses of the important parameters(ix)
Sections7.5 and 7.7 demonstrate large-scale problems(x)
Section 8 gives a conclusion and outlook on future research(xi)
AppendixA provides a more detailed discussion of the calculation of contact forces for the DEM(xii)
AppendixB provides references to the software used in this study and provides scripts for reproduction of results


## 2. The Finite Element Method

The 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–28], kinematic relations describe the possible movement of such elements. The current configuration x can be described with the help of a time-dependent map (noted as t) from the reference configuration X, which can be obtained from the displacement vector u (additionally see Figure 1):(1)xt=X+ut,the bold letters indicate tensors and vectors.Figure 1
Spatial mapping of cable element kinematics between the nodesnk and nl.


### 2.1. Virtual Work

The momentum equation can be described by applying the divergence theorem [22, 27, 29, 30]:(2)divσ+b=ρu¨,withσ being the first Cauchy stress, the body forcesb, 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:(3)δWint−δWext+δWkin=0.The 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:(4)δWint=∫Ωσ:δedΩ,(5)δWint=∫Ω0S:δEdΩ0,(6)δWext=∫Ωb⋅δudΩ+∫Γσσn⋅δudΓσ,(7)δWkin=∫Ωρu¨⋅δudΩ,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.


### 2.2. Cable Element

With regard to the cable element, a cable or a “truss element formulation” having two nodes (nk and nl) considering geometrical nonlinear behaviour is used. The 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 Fint read for each degree of freedom r as follows:(8)Fr,int=A∫0LS+Spre:∂E∂urdX︸L/2⋅dξ,for the constant crosssection A, the reference length L, and a given prestress Spre. The tangent stiffness matrix K can be expressed as(9)Krs=A⋅L2∫−1+1∂S+Spre∂us:∂E∂ur+S+Spre:∂2E∂ur∂usdξ,Krs=A⋅L2∫−1+1D:∂E∂ur:∂E∂us+S+Spre:∂2E∂ur∂usdξ.The 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]:(10)S=∂Ψsv∂E,D=∂2Ψsv∂E∂E.The prestressSpre is a stress, which is initially applied and added to the internal forces of the structure. It is a second-order 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, form-finding procedures need to be applied beforehand. In contrast to standard structural analysis, form-finding solves an inverse engineering problem. This inverse problem is to find a geometry which corresponds to an equilibrium configuration, whose stress state is identical to the given prestress. The force density method is one approach used to realize this [24, 31]. The 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.The element mass matrixM 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:(11)M=A⋅L2⋅ρ0∫−1+1NξTNξdξ.The damping matrixC 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]):(12)C=τM+κK,with τ and κ as the scaling factors.


## 2.1. Virtual Work

The momentum equation can be described by applying the divergence theorem [22, 27, 29, 30]:(2)divσ+b=ρu¨,withσ being the first Cauchy stress, the body forcesb, 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:(3)δWint−δWext+δWkin=0.The 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:(4)δWint=∫Ωσ:δedΩ,(5)δWint=∫Ω0S:δEdΩ0,(6)δWext=∫Ωb⋅δudΩ+∫Γσσn⋅δudΓσ,(7)δWkin=∫Ωρu¨⋅δudΩ,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.


## 2.2. Cable Element

With regard to the cable element, a cable or a “truss element formulation” having two nodes (nk and nl) considering geometrical nonlinear behaviour is used. The 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 Fint read for each degree of freedom r as follows:(8)Fr,int=A∫0LS+Spre:∂E∂urdX︸L/2⋅dξ,for the constant crosssection A, the reference length L, and a given prestress Spre. The tangent stiffness matrix K can be expressed as(9)Krs=A⋅L2∫−1+1∂S+Spre∂us:∂E∂ur+S+Spre:∂2E∂ur∂usdξ,Krs=A⋅L2∫−1+1D:∂E∂ur:∂E∂us+S+Spre:∂2E∂ur∂usdξ.The 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]:(10)S=∂Ψsv∂E,D=∂2Ψsv∂E∂E.The prestressSpre is a stress, which is initially applied and added to the internal forces of the structure. It is a second-order 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, form-finding procedures need to be applied beforehand. In contrast to standard structural analysis, form-finding solves an inverse engineering problem. This inverse problem is to find a geometry which corresponds to an equilibrium configuration, whose stress state is identical to the given prestress. The force density method is one approach used to realize this [24, 31]. The 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.The element mass matrixM 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:(11)M=A⋅L2⋅ρ0∫−1+1NξTNξdξ.The damping matrixC 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]):(12)C=τM+κK,with τ and κ as the scaling factors.


## 3. The Discrete Element Method

This section describes the fundamentals of the discrete element method (DEM), with a focus towards the DEM-FEM coupling (following Sections4–7).The 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. The mass and centre of gravity are described within the cluster shape and independent of the masses of the particles (for further information, see [33]).The basic steps in a DEM simulation are as follows (also [21]):(i)
Contact detection(ii)
Evaluation of contact forces(iii)
Integration of motionThese steps will be briefly described in the following.


### 3.1. Contact Detection

Within the scope of this publication, certain contact scenarios are of interest. The centre of spheres is described withCi, and their respective radius is Ri. The distance to contacted edges is deciphered with di (Figure 2).Two spheresCi and Cj:(13)Ci−Cj<Ri+Rj.SphereCi and vertex nk:(14)Ci−nk<Ri.SphereCi and di to edge:(15)di<Ri.Figure 2
Spherical DEM particlei in contact with cable structure between the nodes nk and nl.


### 3.2. Evaluation of Contact Forces

The 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. The coefficient of restitution (COR), ε⊥ (see equation (A.7)), is an essential part in the contact and is further investigated in Section 7.4.


### 3.3. Integration of Motion

The integration of motion is described by Newton’s second law of motion. The translational accelerationu¨ is described via the force F and the mass of the particle m. The angular acceleration ω˙ is expressed by the moments T and the tensor of rotational inertia I [21]:(16)mu¨=F,Iω˙=T.The forces and torques at each particlei are described as follows [6, 19, 21]:(17)Fi=Fiext+∑j=1nFij︷Fij=Fnnij+Fttij+Fidamp,Ti=Tiext+∑j=1nrcij×Fij+Tidamp.The resulting forces and torques depend on the following components:Fiext,Tiext: external loads (e.g., gravity)Fij: contact interactions between neighbouring particles or boundaries as the contact Fcontact (see Appendix A for description of the components Fn and Ft)Fidamp,Tidamp: external damping boundariesnij,tij: normal and tangential vectors at the respective contact pointrcij: vector connecting the particle centre and the respective contact point


## 3.1. Contact Detection

Within the scope of this publication, certain contact scenarios are of interest. The centre of spheres is described withCi, and their respective radius is Ri. The distance to contacted edges is deciphered with di (Figure 2).Two spheresCi and Cj:(13)Ci−Cj<Ri+Rj.SphereCi and vertex nk:(14)Ci−nk<Ri.SphereCi and di to edge:(15)di<Ri.Figure 2
Spherical DEM particlei in contact with cable structure between the nodes nk and nl.


## 3.2. Evaluation of Contact Forces

The 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. The coefficient of restitution (COR), ε⊥ (see equation (A.7)), is an essential part in the contact and is further investigated in Section 7.4.


## 3.3. Integration of Motion

The integration of motion is described by Newton’s second law of motion. The translational accelerationu¨ is described via the force F and the mass of the particle m. The angular acceleration ω˙ is expressed by the moments T and the tensor of rotational inertia I [21]:(16)mu¨=F,Iω˙=T.The forces and torques at each particlei are described as follows [6, 19, 21]:(17)Fi=Fiext+∑j=1nFij︷Fij=Fnnij+Fttij+Fidamp,Ti=Tiext+∑j=1nrcij×Fij+Tidamp.The resulting forces and torques depend on the following components:Fiext,Tiext: external loads (e.g., gravity)Fij: contact interactions between neighbouring particles or boundaries as the contact Fcontact (see Appendix A for description of the components Fn and Ft)Fidamp,Tidamp: external damping boundariesnij,tij: normal and tangential vectors at the respective contact pointrcij: vector connecting the particle centre and the respective contact point


## 4. 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.


### 4.1. 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:(18)FcontactuΩS,Γt,u˙ΩS,Γt,uPt,u˙Pt−FintΩS,ΓuΩS,Γt,u˙ΩS,Γt,u¨ΩS,Γt=0,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 Fcontact of the particles, which are dependent on its displacements and the velocities, need to be in equilibrium with the internal forces FintΩS,Γ of the structure, which are dependent on its displacements, velocities, and accelerations (equation (8) and additional including inertia effects). The equilibrium of both DEM and FEM simulation is graphically depicted in Figure 3.Figure 3
Domain definition and force equilibrium in DEM and FEM simulation.The 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.


### 4.2. 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. The 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. This procedure is outlined in Figure 4.Figure 4
Transfer of forces and displacements and velocities between different applications.Following this, the contact forcesFcontact 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:(19)FextΩS,Γ=FcontactuΩD,Γt,u˙ΩD,Γt,uPt,u˙Pt.The equilibrium within the structural mechanics problem is given as follows:(20)FextΩS,Γ−FintΩS,ΓuΩS,Γt,u˙ΩS,Γt,u¨ΩS,Γt=0.After solving both domains, the two interface conditions, for the displacements and the velocities between both fields, are not fulfilled anymore:(21)uΩD,Γt−uΩS,Γt=0,(22)u˙ΩD,Γt−u˙ΩS,Γt=0.Resulting from this, the contact forcesFcontact 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:(23)FcontactuΩD,Γt,u˙ΩD,Γt,uPt,u˙Pt−FcontactuΩS,Γt,u˙ΩS,Γt,uPt,u˙Pt=0.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 and makes the solution unstable. To solve this problem, a possible approach is presented inSection 6.


### 4.3. 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. This value must be defined by the user and heavily influences the stability of the coupled simulation (see example in Section 7.4).The 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.


### 4.4. 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 cable-like structures, the DEM wall conditionΩD discretization and the FEM ΩS discretization on Γ must exactly coincide (conforming meshes). The 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.Figure 5
Modelling of impact with cable-structured boundaries. (a) DiscretizedΩD. (b) Discretized ΩS.
(a)(b)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.


## 4.1. 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:(18)FcontactuΩS,Γt,u˙ΩS,Γt,uPt,u˙Pt−FintΩS,ΓuΩS,Γt,u˙ΩS,Γt,u¨ΩS,Γt=0,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 Fcontact of the particles, which are dependent on its displacements and the velocities, need to be in equilibrium with the internal forces FintΩS,Γ of the structure, which are dependent on its displacements, velocities, and accelerations (equation (8) and additional including inertia effects). The equilibrium of both DEM and FEM simulation is graphically depicted in Figure 3.Figure 3
Domain definition and force equilibrium in DEM and FEM simulation.The 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.


## 4.2. 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. The 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. This procedure is outlined in Figure 4.Figure 4
Transfer of forces and displacements and velocities between different applications.Following this, the contact forcesFcontact 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:(19)FextΩS,Γ=FcontactuΩD,Γt,u˙ΩD,Γt,uPt,u˙Pt.The equilibrium within the structural mechanics problem is given as follows:(20)FextΩS,Γ−FintΩS,ΓuΩS,Γt,u˙ΩS,Γt,u¨ΩS,Γt=0.After solving both domains, the two interface conditions, for the displacements and the velocities between both fields, are not fulfilled anymore:(21)uΩD,Γt−uΩS,Γt=0,(22)u˙ΩD,Γt−u˙ΩS,Γt=0.Resulting from this, the contact forcesFcontact 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:(23)FcontactuΩD,Γt,u˙ΩD,Γt,uPt,u˙Pt−FcontactuΩS,Γt,u˙ΩS,Γt,uPt,u˙Pt=0.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 and makes the solution unstable. To solve this problem, a possible approach is presented inSection 6.


## 4.3. 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. This value must be defined by the user and heavily influences the stability of the coupled simulation (see example in Section 7.4).The 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.


## 4.4. 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 cable-like structures, the DEM wall conditionΩD discretization and the FEM ΩS discretization on Γ must exactly coincide (conforming meshes). The 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.Figure 5
Modelling of impact with cable-structured boundaries. (a) DiscretizedΩD. (b) Discretized ΩS.
(a)(b)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.


## 5. Staggered Weak Coupling

The fundamental idea of the weak coupling (sometimes also called explicit coupling [24]) follows a single exchange of coupling data in each time step. The communication pattern is depicted in Figure 6. The important steps at each time, including this communication pattern, can be summarized as follows:(1)
Solve DEM (results:uP,u˙P,andFcontact)(2)
MapFcontact from DEM to Structure(3)
Solve Structure (results:uΩS,Γ,u˙ΩS,Γ,andu¨ΩS,Γ)(4)
Mapdisplacements and velocities from Structure to DEM(5)
Advance in time (not explicitly shown)Figure 6
Staggered weak coupling procedure, between DEM and FEM.The interface variables are accordingly updated (see Steps (2) and (4)):Displacement:(24)uk+1ΩD,Γ=ukΩS,Γ.Velocity:(25)u˙k+1ΩD,Γ=u˙kΩS,Γ.Contact force:(26)Fk+1extΩS,Γ=Fkcontact..This 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. Therefore, 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 Sections7.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. The 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 Algorithm1. In order to facilitate the reproduction of the results, the Python script used, including all comments, is provided in Appendix B.Algorithm 1:Weak coupling.
(1)
Initialize(2)
Whiletime < end_time do(3)
Ifparticle_near_wallthen(4)
use predefined time step(5)
Else(6)
use increased time step(7)
t=t+Δt(8)
Search nearest neighbours and find contact ⟶ equations (13)–(15)(9)
Calculate contact forces ⟶ equation (A.5)(10)
Time integration of DEM part ⟶ equation (16)(11)
Mapforces on Γ from ΩD to ΩS(12)
Ifforces≠0.0then(13)
Solve structure (FEM)(14)
Mapvelocity and displacement on Γ from ΩS to ΩD(15)
Update position ofΩD to ΩS ⟶ equation (1)(16)
FinalizeIn this procedure, two additional features will be discussed. They are independent of the coupling approach used but improve the performance significantly. They are added in Algorithm1 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–15): this is an additional check used only to solve the structure when contact forces are present. This 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.


## 6. 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. This 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. This requires a Gauss–Seidel loop between DEM and FEM, which might need to be solved multiple times within one time step [23, 24, 38, 39]. This 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 ε. The residual formulation is defined by using equation (27)The steps of this approach are shown in Figure7 and summarized in the following, using the respective numbering in the abovementioned Figure 7:(1)
Solve DEM (results:uP,u˙P,andFcontact)(2)
MapFcontact from DEM to Structure(3)
Solve Structure (results:uΩS,Γ,u˙ΩS,Γ,andu¨ΩS,Γ)(4)
Mapdisplacements 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 timeFigure 7
Strong coupling communication diagram.The weak coupling algorithm, described in the preceding Section5 expresses single iteration in the strong coupling scheme (Steps (1)–(4)). The 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. This is especially due to the property that many coupling iterations are typically not required throughout the simulation, but only at specific time steps. The comparison of the two procedures, including a view on the performance, is outlined in Section 7.The residual criteria within the strong coupling loop are defined by(27)ϵ<rkneq,where ϵ is the user-defined breaking tolerance. It is checked after each iteration k by scaling the norm of the residuum rk with the square root of the number of degrees of freedom neq 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. The 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:Displacement residuum:(28)rku=ukΩS,Γ−uk−1ΩS,Γ.Velocity residuum:(29)rku˙=u˙kΩS,Γ−u˙k−1ΩS,Γ.Contact force residuum:(30)rkF=Fkcontact−Fk−1contact.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. The outcome is that this permits a faster interface convergence. The so-calledconvergence 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. The relaxation is done w.r.t. the residual (equations (28)–(30)), respectively:Relaxed displacements:(31)ukrelΩS,Γ=uk−1relΩS,Γ+αkurku.Relaxed velocities:(32)u˙krelΩS,Γ=u˙k−1relΩS,Γ+αku˙rku˙.Relaxed contact forces:(33)Fkcontact,rel=Fk−1contact,rel+αkFrkF.Each variable is subsequently updated from the previous solution (Stepk−1) using the respective interface residuum scaled by relaxation factor α.There are different approaches to obtain the scaling factorα [39]. The relaxation factor can be set to a user-defined constant value, which is very simple and helps to improve the quality of the simulation. Another approach is to use the Aitken method [41]. The Aitken method optimizes α in every iteration w.r.t. the current residuum rk and the previous residuum rk−1:(34)αkr=−αk−1rk−1⋅rk−rk−1rk−rk−1,respectively, αu=αru,αu˙=αru˙,andαF=αrF. The 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. Thus, the resulting residuum for both relaxing procedures is given as follows:Displacement and velocity residua:(35)rk=maxrku,rku˙.Contact force residuum:(36)rk=rkF.The interface variables are updated accordingly (see Steps 2 and 4 in Figure7). The following variables are exchanged within the interface:Without relaxationDisplacements:(37)uk+1ΩD,Γ=ukΩS,Γ.Velocities:(38)u˙k+1ΩD,Γ=u˙kΩS,Γ.Contact forces:(39)Fk+1extΩS,Γ=Fkcontact.With relaxationDisplacements:(40)uk+1ΩD,Γ=ukrelΩS,Γ.Velocities:(41)u˙k+1ΩD,Γ=u˙krelΩS,Γ.Contact forces:(42)Fk+1extΩS,Γ=Fkcontact,rel.In summary, both solution strategies are described within Algorithms2 and 3 in pseudocode. They are both further elaborated on in Appendix B.3.Algorithm 2:Strong coupling: relaxed displacement and velocity.
(1)
Initialize(2)
Whiletime < end_timedo(3)
Whileinterface_res ≥ tolerance_interfacedo(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)
Mapforces on Γ from ΩD to ΩS(8)
Solve structure (FEM)(9)
Mapvelocity and displacement on Γ from ΩS to ΩD(10)
Calculate residuals forvelocity and displacement ⟶ equations (28) and (29)(11)
Relaxvelocity and displacement ⟶ equations (31) and (32)(12)
Update position ofΩD ⟶ equation (1)(13)
interface_res = max(displacement_residual, velocity_residual)(14)
Update position ofΩS ⟶ equation (1)(15)
FinalizeAlgorithm 3:Strong coupling: relaxed forces.
(1)
Initialize(2)
Whiletime < end_timedo(3)
Whileinterface_res ≥ tolerance_interfacedo(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)
Mapforces on Γ from ΩDto ΩS(8)
Calculate residuals forforces ⟶ equation (30)(9)
Relaxforces ⟶ equation (33)(10)
Solve structure (FEM)(11)
Mapvelocity and displacement on Γ from Ωsto ΩD(12)
Update position ofΩD ⟶ equation (1)(13)
interface_res = force_residual(14)
Update position ofΩS ⟶ equation (1)(15)
Finalize


## 7. Systematic Assessment of the DEM-FEM Coupling

This 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. The examples show problems of impacting objects on highly flexible lightweight cable structures, such as protection nets. These interaction problems typically have numerical stability issues within the simulations, as the net structures have a low mass, whereas the rocks are typically heavy. This instability leads to the problem that especially when the first impact occurs, the forces might become very large. Thus, 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 problem7.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). The proper choice of a relaxation factor is further discussed in the example in 7.3. The 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.


### 7.1. 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. The setup of this academic example can be found in Figure8(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.Figure 8
Problem setup, including all geometrical and physical input parameters for the simulation. Results after impact with time step ofΔt=10−2s. (a) Problem setup. (b) Weak coupling after impact. (c) Strong coupling after impact.
(a)(b)(c)Within empirical tests, the time stepΔt=10−3s is found to be the highest possible time step for which the weak coupling algorithm can resolve to an appropriate solution. 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−2s. 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. This leads to a nonphysical behaviour of the coupled problem, as shown in Figure 9.Figure 9
Comparison of weak coupling approach withΔt=10−2s and Δt=10−3s against the strong coupling approach with Δt=10−2s for the vertical displacements of node (A) (Figure 8(a)).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−2s.In the following, the time step of the first contact is discussed in detail. It can be seen (Figure10(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). This exemplarily demonstrates the advantages of the strong coupling scheme, presented in this article. The large discrepancy in the contact force would lead to an unstable coupled simulation when using a standard weak coupling algorithm. The 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)).Figure 10
Simulation results within one time step over inner coupling iterations. (a) Interface displacement. (b) Interface velocity. (c) Displacement of impacting sphere. (d) Contact forces onΩD.
(a)(b)(c)(d)As an example, the converging contact force for each iteration in time stept=1.2s 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−3s, which is the limit for which the weak coupling approach still converges.Figure 11
Interface simulation results. (a) Contact force att=1.2s. (b) Number of interface iterations per time step.
(a)(b)


### 7.2. 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 (Spre=106N/m2) and a homogeneous Young's Modulus (E=109N/m2) in the cable structure is used with the following changes for the impacting sphere: R=0.12m andρ=3.5⋅104kg/m3. The result of the transient analysis will be compared to the static solution, considering the particle as an external force. This static force is defined as follows:(43)Fstatic=43πRsphere3︸Vsphere⋅ρsphere︷msphere⋅9.81ms2︸gravity.The resulting static deflection of Point A (Figure8(a)) is shown in Figure 12. This comparison proves that the transient analysis approaches the static solution after a certain time.Figure 12
Comparison of weak coupling approach withΔt=10−2s and Δt=10−3s against the strong coupling approach with Δt=10−2s and Δt=3⋅10−2s 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.Furthermore, the sensitivity of the time step within each coupling algorithm is also studied in this example. The results of all solutions are presented and compared in Figure12. It shows that the weak coupling approach provides an accurate performance for a time step of Δt=10−3s, whereas the solution for Δt=10−2s is very unstable. It turns out that for large time steps, the result oscillates around the expected solutions.The measured solutions for time steps ofΔt=10−2s and Δt=3⋅10−2s 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.Figure 13
Number of interface iterations per time step, for either (a) relaxed displacements and velocities or (b) relaxed forces.
(a)(b)The influence if either displacements and velocities or forces are relaxed is examined in the following. Both options are described in Algorithms2 and 3, respectively.Comparing Figures13(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 followingSection 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.


### 7.3. 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. The same problem setup as in Section 7.2 is used.As Figure14 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 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.Figure 14
Comparison of displacement over time for different relaxation factorsα. 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.Table 1
Comparison of computation time with different relaxations and theAitken relaxation.
Relaxation factorαRelative computation time (%)0.1100.00.254.10.520.40.728.5Aitken18.7


### 7.4. 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–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.Figure 15
Damping coefficientζε⊥, see [34, 35] for more information. (a) Influence of COR w.r.t. equation (A.7). (b) Influence of COR w.r.t. equation (A.8).
(a)(b)As can be seen in Figure16(a) (zoomed in Figure 16(b)), the interface coupling becomes unstable as soon as the COR reaches a small value. This 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 17 visualizes the progression of the maximum interface iterations over the simulation time and indicates the advantages of the proposed coupling algorithm. The large number of iterations at the time of first contact t≈0.25s calls for a small time step due to the increased difference in interface velocities. This can be overcome with the help of 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.Figure 16
Comparison of different COR values. (a)t0=0.0−t1=2.0. (b) Scaling to t0=0.4−t1=1.4.
(a)(b)Figure 17
Strong coupling interface iterations forε⊥=0.2.


### 7.5. Practical Application: Angled Protection Net

One prominent practical application case of highly flexible structures can be found in mountainous regions. As an 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 Figure18(a).Figure 18
Picture of protection net and corresponding numerical model, including the net as cable structure and two impacting rocks modelled as single DEM particles. The one particle is chosen to be so small that it falls through the net, to verify the correct coupling detection. (a) Route Chalais-Vercorin, Valais [45]. (b) Numerical DEM-FEM coupled model.
(a)(b)To test the limits of the presented algorithms, in this study, the same system as shown in Figure18(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.Table 2
Properties of FEM and DEM parts of example inSection 7.5.
FEMDEMEN/m27⋅105EN/m21⋅105Am21.26⋅10−5u˙0m/s0.0,−5.54,−5.54ρkg/m37.85⋅103ρkg/m32.5⋅103Spre0Rim0.21,0.05τ20ν2⋅10−1κ0ε⊥2⋅10−1Using a time step of2⋅10−4s, the different behaviours after impact are presented in the incidental Figures 19(a) and 19(b).Figure 19
Comparison of the (a) weak coupling to the (b) strong coupling approach with the same time stepΔt=2⋅10−4s.
(a)(b)Similar to the example fromSection 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)).Figure 20
Simulation results. (a) Weak vs. strong coupling-displacement of centre node. (b) Interface iterations over simulation time.
(a)(b)The considerably large number of interface iterations (Figure20(b)) is a result of the system setup. This example 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.


### 7.6. Arbitrary Boundary Conditions

Another advantage of this procedure is the possibility to use arbitrary boundary conditions for the problem. As can be seen in Figures21(a)–22(e), the arbitrary triangular and quadrilateral meshes can be used to simulate any boundary condition, for instance, a mountainous region. As it is the state of the art in industrial applications, ∗.stl files can be used. If only point clouds are available, standard tools can be used to create a triangulated mesh. The structural part can be subsequently put into this domain to easily capture the interaction of different terrain models and loading scenarios.Figure 21
Plane boundary rigid with one impacting sphere and two protection nets. The first net (right) is not connected on the lower edge, and the other is coupled to the ground. (a)t=1.3s. (b) t=1.5s. (c) t=1.9s. (d) t=2.2s. (e) t=2.4s. (f) t=2.5s. (g) t=3.6s. (h) t=4.4s. (i) t=4.7s.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 22
Triangulated terrain model, showing the downward rolling of the rocks, which are being slowed and stopped by the protection net. (a)t=1.2s. (b) t=4.5s. (c) t=6.8s. (d) Net impact t=6.8s. (e) Net impact: side view t=6.8s.
(a)(b)(c)(d)(e)As an example, Figures21(a)–21(i) demonstrate a plane boundary with a curtain-like structure in the middle and a cable net protection net at the end of the slope.In contrast to the plane boundary, Figures22(a)–22(e) show the use of an arbitrarily shaped boundary, obtained from a ∗.stl mesh. The effortless integration of a deformable FEM structure into the arbitrary boundary is indicated in this example.


### 7.7. 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-of-the-art rockfall protection simulations as discussed in [1, 47].Figure 23
Sliding node within the FEM model. It allows us to model a more accurate behaviour of the cable structure nets. (a)t=0.0s. (b) t=2.7s. (c) t=4.5s.
(a)(b)(c)Figure 24
Impact problem with more complex shape of impacting objects. Moving impact rocks are discretized as particles which are clustered to allow any arbitrary shapes. (a)t=0.0s. (b) t=1.3s. (c) t=3.0s. (d) Impact of DEM clusters, t=1.3s.
(a)(b)(c)(d)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. The 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.The 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.Figure 25
Small particle is penetrating and falling through the coarse cable net. Big particle collides with cable structure and is stopped and thrown back by the protection structure. (a)t=0.0s. (b) t=0.05s. (c) t=0.1s. (d) t=0.15s. (e) t=0.2s.Figure 26
Protection net test setup in Walenstadt, Switzerland. (a) Testing area in Walenstadt, St. Gallen, Switzerland [45]. (b) FEM mesh.
(a)(b)


## 7.1. 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. The setup of this academic example can be found in Figure8(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.Figure 8
Problem setup, including all geometrical and physical input parameters for the simulation. Results after impact with time step ofΔt=10−2s. (a) Problem setup. (b) Weak coupling after impact. (c) Strong coupling after impact.
(a)(b)(c)Within empirical tests, the time stepΔt=10−3s is found to be the highest possible time step for which the weak coupling algorithm can resolve to an appropriate solution. 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−2s. 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. This leads to a nonphysical behaviour of the coupled problem, as shown in Figure 9.Figure 9
Comparison of weak coupling approach withΔt=10−2s and Δt=10−3s against the strong coupling approach with Δt=10−2s for the vertical displacements of node (A) (Figure 8(a)).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−2s.In the following, the time step of the first contact is discussed in detail. It can be seen (Figure10(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). This exemplarily demonstrates the advantages of the strong coupling scheme, presented in this article. The large discrepancy in the contact force would lead to an unstable coupled simulation when using a standard weak coupling algorithm. The 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)).Figure 10
Simulation results within one time step over inner coupling iterations. (a) Interface displacement. (b) Interface velocity. (c) Displacement of impacting sphere. (d) Contact forces onΩD.
(a)(b)(c)(d)As an example, the converging contact force for each iteration in time stept=1.2s 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−3s, which is the limit for which the weak coupling approach still converges.Figure 11
Interface simulation results. (a) Contact force att=1.2s. (b) Number of interface iterations per time step.
(a)(b)


## 7.2. 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 (Spre=106N/m2) and a homogeneous Young's Modulus (E=109N/m2) in the cable structure is used with the following changes for the impacting sphere: R=0.12m andρ=3.5⋅104kg/m3. The result of the transient analysis will be compared to the static solution, considering the particle as an external force. This static force is defined as follows:(43)Fstatic=43πRsphere3︸Vsphere⋅ρsphere︷msphere⋅9.81ms2︸gravity.The resulting static deflection of Point A (Figure8(a)) is shown in Figure 12. This comparison proves that the transient analysis approaches the static solution after a certain time.Figure 12
Comparison of weak coupling approach withΔt=10−2s and Δt=10−3s against the strong coupling approach with Δt=10−2s and Δt=3⋅10−2s 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.Furthermore, the sensitivity of the time step within each coupling algorithm is also studied in this example. The results of all solutions are presented and compared in Figure12. It shows that the weak coupling approach provides an accurate performance for a time step of Δt=10−3s, whereas the solution for Δt=10−2s is very unstable. It turns out that for large time steps, the result oscillates around the expected solutions.The measured solutions for time steps ofΔt=10−2s and Δt=3⋅10−2s 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.Figure 13
Number of interface iterations per time step, for either (a) relaxed displacements and velocities or (b) relaxed forces.
(a)(b)The influence if either displacements and velocities or forces are relaxed is examined in the following. Both options are described in Algorithms2 and 3, respectively.Comparing Figures13(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 followingSection 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.


## 7.3. 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. The same problem setup as in Section 7.2 is used.As Figure14 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 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.Figure 14
Comparison of displacement over time for different relaxation factorsα. 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.Table 1
Comparison of computation time with different relaxations and theAitken relaxation.
Relaxation factorαRelative computation time (%)0.1100.00.254.10.520.40.728.5Aitken18.7


## 7.4. 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–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.Figure 15
Damping coefficientζε⊥, see [34, 35] for more information. (a) Influence of COR w.r.t. equation (A.7). (b) Influence of COR w.r.t. equation (A.8).
(a)(b)As can be seen in Figure16(a) (zoomed in Figure 16(b)), the interface coupling becomes unstable as soon as the COR reaches a small value. This 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 17 visualizes the progression of the maximum interface iterations over the simulation time and indicates the advantages of the proposed coupling algorithm. The large number of iterations at the time of first contact t≈0.25s calls for a small time step due to the increased difference in interface velocities. This can be overcome with the help of 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.Figure 16
Comparison of different COR values. (a)t0=0.0−t1=2.0. (b) Scaling to t0=0.4−t1=1.4.
(a)(b)Figure 17
Strong coupling interface iterations forε⊥=0.2.


## 7.5. Practical Application: Angled Protection Net

One prominent practical application case of highly flexible structures can be found in mountainous regions. As an 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 Figure18(a).Figure 18
Picture of protection net and corresponding numerical model, including the net as cable structure and two impacting rocks modelled as single DEM particles. The one particle is chosen to be so small that it falls through the net, to verify the correct coupling detection. (a) Route Chalais-Vercorin, Valais [45]. (b) Numerical DEM-FEM coupled model.
(a)(b)To test the limits of the presented algorithms, in this study, the same system as shown in Figure18(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.Table 2
Properties of FEM and DEM parts of example inSection 7.5.
FEMDEMEN/m27⋅105EN/m21⋅105Am21.26⋅10−5u˙0m/s0.0,−5.54,−5.54ρkg/m37.85⋅103ρkg/m32.5⋅103Spre0Rim0.21,0.05τ20ν2⋅10−1κ0ε⊥2⋅10−1Using a time step of2⋅10−4s, the different behaviours after impact are presented in the incidental Figures 19(a) and 19(b).Figure 19
Comparison of the (a) weak coupling to the (b) strong coupling approach with the same time stepΔt=2⋅10−4s.
(a)(b)Similar to the example fromSection 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)).Figure 20
Simulation results. (a) Weak vs. strong coupling-displacement of centre node. (b) Interface iterations over simulation time.
(a)(b)The considerably large number of interface iterations (Figure20(b)) is a result of the system setup. This example 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.


## 7.6. Arbitrary Boundary Conditions

Another advantage of this procedure is the possibility to use arbitrary boundary conditions for the problem. As can be seen in Figures21(a)–22(e), the arbitrary triangular and quadrilateral meshes can be used to simulate any boundary condition, for instance, a mountainous region. As it is the state of the art in industrial applications, ∗.stl files can be used. If only point clouds are available, standard tools can be used to create a triangulated mesh. The structural part can be subsequently put into this domain to easily capture the interaction of different terrain models and loading scenarios.Figure 21
Plane boundary rigid with one impacting sphere and two protection nets. The first net (right) is not connected on the lower edge, and the other is coupled to the ground. (a)t=1.3s. (b) t=1.5s. (c) t=1.9s. (d) t=2.2s. (e) t=2.4s. (f) t=2.5s. (g) t=3.6s. (h) t=4.4s. (i) t=4.7s.
(a)(b)(c)(d)(e)(f)(g)(h)(i)Figure 22
Triangulated terrain model, showing the downward rolling of the rocks, which are being slowed and stopped by the protection net. (a)t=1.2s. (b) t=4.5s. (c) t=6.8s. (d) Net impact t=6.8s. (e) Net impact: side view t=6.8s.
(a)(b)(c)(d)(e)As an example, Figures21(a)–21(i) demonstrate a plane boundary with a curtain-like structure in the middle and a cable net protection net at the end of the slope.In contrast to the plane boundary, Figures22(a)–22(e) show the use of an arbitrarily shaped boundary, obtained from a ∗.stl mesh. The effortless integration of a deformable FEM structure into the arbitrary boundary is indicated in this example.


## 7.7. 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-of-the-art rockfall protection simulations as discussed in [1, 47].Figure 23
Sliding node within the FEM model. It allows us to model a more accurate behaviour of the cable structure nets. (a)t=0.0s. (b) t=2.7s. (c) t=4.5s.
(a)(b)(c)Figure 24
Impact problem with more complex shape of impacting objects. Moving impact rocks are discretized as particles which are clustered to allow any arbitrary shapes. (a)t=0.0s. (b) t=1.3s. (c) t=3.0s. (d) Impact of DEM clusters, t=1.3s.
(a)(b)(c)(d)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. The 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.The 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.Figure 25
Small particle is penetrating and falling through the coarse cable net. Big particle collides with cable structure and is stopped and thrown back by the protection structure. (a)t=0.0s. (b) t=0.05s. (c) t=0.1s. (d) t=0.15s. (e) t=0.2s.Figure 26
Protection net test setup in Walenstadt, Switzerland. (a) Testing area in Walenstadt, St. Gallen, Switzerland [45]. (b) FEM mesh.
(a)(b)


## 8. Conclusions and Outlook

The 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.The 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. Thus, to overcome this problem, an additional iteration between the physics was explained in Section 6. This 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.The novel approaches make it possible to efficiently simulate the correct behaviour of complex existing structures. The example in7.5 shows net structures interacted with rocks which are based on existing structures in the Austrian and Swiss Alps. The 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 inSection 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)). This 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.



---

*Source: 5135194-2020-11-21.xml*