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.
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)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
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 can be described with the help of a time-dependent map (noted as ) from the reference configuration , which can be obtained from the displacement vector (additionally see Figure 1):the bold letters indicate tensors and vectors.
2.1. Virtual Work
The momentum equation can be described by applying the divergence theorem [22, 27, 29, 30]:with being the first Cauchy stress, the body forces, the acceleration , and the density . Multiplying the momentum equation (2) with the virtual displacement field and integrating over the physical domain leads to the weak form, respectively:
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:with the Almansi strain , the normal direction of , , 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 and the Green–Lagrange strain tensor within the initial domain .
2.2. Cable Element
With regard to the cable element, a cable or a “truss element formulation” having two nodes ( and ) 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 is used to describe the geometry [25] with linear shape functions.
By virtue of equation (5), the internal forces read for each degree of freedom as follows:for the constant crosssection , the reference length , and a given prestress . The tangent stiffness matrix can be expressed as
The constitutive law used in this work is the St. Venant–Kirchhoff hyperelastic material model [29, 30] described by the strain-energy functional . Its first and second derivatives w.r.t. describe the PK2-stress and the material tangent modulus , respectively [30]:
The prestress 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 matrix for a dynamic simulation results from the virtual kinetic work as depicted in equation (7). With being the shape functions and the reference density , is defined as follows:
The damping matrix is expressed via the Rayleigh damping as a linear combination of the mass matrix and the stiffness matrix (for more information, see [32]):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 Sections 4–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 motion
These 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 with , and their respective radius is . The distance to contacted edges is deciphered with (Figure 2). Two spheres and : Sphere and vertex : Sphere and to edge:
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 acceleration is described via the force and the mass of the particle . The angular acceleration is expressed by the moments and the tensor of rotational inertia [21]:
The forces and torques at each particle are described as follows [6, 19, 21]:
The resulting forces and torques depend on the following components: : external loads (e.g., gravity) : contact interactions between neighbouring particles or boundaries as the contact (see Appendix A for description of the components and ) : external damping boundaries : normal and tangential vectors at the respective contact point : 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:where describes the displacements, represents the velocities, represents the structural domain, while of includes all nodes which are in interaction with the DEM particles, .
With respect to equation (18), the contact forces of the particles, which are dependent on its displacements and the velocities, need to be in equilibrium with the internal forces 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.
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 . The wall is used to calculate contact forces with the DEM particles , 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 . This procedure is outlined in Figure 4.
Following this, the contact forces are now dependent on the displacements and velocities of and not directly on and are defined as the external forces coming from the DEM:
The 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 computed within and the contact forces which would be computed within 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 and makes the solution unstable. To solve this problem, a possible approach is presented in Section 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 discretization and the FEM 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.
(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: )(2)Map from DEM to Structure(3)Solve Structure (results: )(4)Map displacements and velocities from Structure to DEM(5)Advance in time (not explicitly shown)
The interface variables are accordingly updated (see Steps (2) and (4)): Displacement: Velocity: Contact force:
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 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. 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 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. They are independent of the coupling approach used but improve the performance significantly. They 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 (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 Figure 7 and summarized in the following, using the respective numbering in the abovementioned Figure 7:(1)Solve DEM (results: )(2)Map from DEM to Structure(3)Solve Structure (results: )(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
The weak coupling algorithm, described in the preceding Section 5 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 bywhere is the user-defined breaking tolerance. It is checked after each iteration by scaling the norm of the residuum with the square root of the number of degrees of freedom 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 , the residuum of each variable can be noted as follows: Displacement residuum: 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. The outcome is that this permits a faster interface convergence. The 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. The relaxation is done w.r.t. the residual (equations (28)–(30)), respectively: Relaxed displacements: Relaxed velocities: Relaxed contact forces:
Each variable is subsequently updated from the previous solution (Step ) 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 and the previous residuum :respectively, . 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: Contact force residuum:
The interface variables are updated accordingly (see Steps 2 and 4 in Figure 7). The 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. They are both further elaborated on in Appendix B.3.
|
|
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 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). 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 Figure 8(a), with 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.
(a)
(b)
(c)
Within empirical tests, the time step 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 . 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 . 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.
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 .
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). 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 and do converge to a common value to fulfil the interface displacement/velocity equilibrium (equations (21) and (22)).
(a)
(b)
(c)
(d)
As an example, the converging contact force for each iteration in time step 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 , which is the limit for which the weak coupling approach still converges.
(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 () and a homogeneous Young's Modulus () in the cable structure is used with the following changes for the impacting sphere: . 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:
The resulting static deflection of Point A (Figure 8(a)) is shown in Figure 12. This 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. The 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 , whereas the solution for 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 and 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.
(a)
(b)
The 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.
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 (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 Figure 14 shows, a constant relaxation factor can be used as long as it is smaller than . 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.
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.
(a)
(b)
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. 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 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.
(a)
(b)
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 Figure 18(a).
(a)
(b)
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 and a high impact velocity are chosen to introduce even more difficulties due to an increased contact force.
Using a time step of , the different behaviours after impact are presented in the incidental Figures 19(a) and 19(b).
(a)
(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)).
(a)
(b)
The considerably large number of interface iterations (Figure 20(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 Figures 21(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.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(a)
(b)
(c)
(d)
(e)
As an example, Figures 21(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, Figures 22(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].
(a)
(b)
(c)
(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 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.
(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 in 7.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 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)). 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.
Appendix
A. DEM Force Derivation
A detailed description of the evaluation of forces described in Section 3 is provided in order to further discuss the necessary quantities in the underlying coupling scheme. As soon as a contact is detected, the forces can be evaluated using various contact laws and rheological models in which the normal indentation [21] is as follows:and its time derivative , the normal vector , and the increment of tangential displacement are used.
In order to obtain , the tangent unit vector must be derived by splitting the velocity at the contact point into normal and tangential components. First, the contact point velocity is expressed with the aid of the respective element velocity , element angular velocity , and the vector connecting the particle centre and the contact point:and then, split upwhich allows us to express using the tangent unit vector, the element displacement , and the element rotation :
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 and tangential contact forces [21, 34] for the case of two spheres colliding,are calculated with the aid of the normal and tangential stiffness and , respectively, and the damping coefficients and [21], considering the maximum tangential force restricted to the Coulomb’s friction limit [34] with the coefficient of friction . The tangential forces are consequently updated from the last step, indicated by the superscript .
(a)
(b)
The material parameters in Table 3 (Young’s modulus , particle mass , shear modulus , and Poisson’s ratio ) are typically obtained by calibration from experiments [49].
As a scaling factor in Table 3, the dashpot coefficient [21, 34] is frequently expressed using the normal coefficient of restitution (COR) as shown in the following equation:
Based on our experience, which is as well as supported by [34, 35], a more realistic modelling of the dashpot coefficient is given by the following equation [34, 35]:
The difference can be observed in Figures 15(a) and 15(b). For further information about this topic, the reader is redirected to [21, 34, 35].
With respect to [21, 35], the COR expresses the ratio between the velocity after and the velocity before impact. A maximum of will model a perfectly elastic impact, whereas 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 [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. The 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
The script to define the problem setup is shown in the following. This 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. This 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. These 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
The 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).
B.3.1. Relax Displacements and Velocities algorithm 7
|
B.3.2. Relax Forces
|
Abbreviations
| FEM: | Finite element method |
| DEM: | Discrete element method |
| COR: | Coefficient of restitution |
| MPM: | Material point method |
| FSI: | Fluid-structure interaction. |
Data Availability
No data were used to support the findings of the study.
Conflicts of Interest
The 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.
Acknowledgments
This work was supported by the Technical University of Munich (TUM).