A unified and modular coupling of particle methods with fem for civil engineering problems

In this work, a modular coupling approach for particle methods with the FEM (finite element method) is presented. The proposed coupled strategy takes advantage from the ability of particle methods of dealing with large displacements and deformations, especially when solving complex fluid–structure and solid–structure interaction problems. The coupling between the FEM and particle methods is done using a co-simulation approach implemented in the open-source Kratos Multiphysics framework. The particle methods considered in this work are the DEM (discrete element method) and the PFEM (particle finite element method). The Lagrangian description of the PFEM is well suited for modeling fluids undergoing large deformations and free-surface motions, and the DEM can be used to simulate rocks, debris and other solid objects. To accelerate the convergence of the coupled strategy, a block Gauss–Seidel algorithm with Aitken relaxation is used. Several numerical examples, with an emphasis on natural hazards, are presented to test and validate the proposed coupled method.


Introduction
Natural hazards modeling in civil engineering is a research area that has drawn increasing attention in the recent years. Due to the climate change trends, extreme natural events, such as floods, landslides, rockfalls, and debris flows, are increasing in frequency every year. For this reason, it is urgent to further investigate modeling approaches to improve the current predictive methodologies of these challenging mul- tiphysics engineering problems. In this problem spectrum, this work analyzes the impact phenomena of big water waves and destructive landslides, such as rockslides or debris flows, against retaining structures. In particular, special focus has been put on protection nets that are used to dissipate the energy from rockfalls, and on membranes and walls subjected to water impacts.
The distinct physics involved in these coupled analyses demands an appropriate computational method selection. Particle methods represent a useful choice for modeling the mentioned natural hazards because they can deal with the large displacements and shape changes undergone by the mobilized material. For example, the discrete element method (DEM) [15], thanks to its efficient contact detection algorithms, is suitable for the simulation of rock assemblies and granular materials. Examples of application of the DEM to the simulation of rockfalls and debris flows can be found in [4,61] and [36,37], respectively. Analogously, free-surface fluid flows, such as tsunami waves or landslides, can be modeled with continuous particle-based methods, such as the material point method (MPM) [27,59,60], the smoothed particle hydrodynamics (SPH) method [13,43,44], or the particle finite element method (PFEM) [23,35,46]. Among the mentioned methods, the PFEM is here employed due to its easier coupling with the finite element method used for the solid structures. The PFEM is a Lagrangian FEM-based method that is capable to track the highly deforming fluid body while avoiding the distortion of the FE mesh thanks to an efficient remeshing technique. Remarkably, the PFEM remeshing allows recognizing not only the evolving fluid free surface but also the contours of the surrounding solid bodies, making the PFEM a well-suited method for fluid-structure interaction (FSI) problems.
For both the protection nets and the flexible membranes, a standard finite element method (FEM) [26,69] is chosen.
The coupling of the FEM with the particle methods is performed in a partitioned way by using an enhanced cosimulation approach, based on the previous works of [8] and [24]. The co-simulation approach has an optimized structure that gives high flexibility when solving coupled problems. The different methods are made to interact with each other through a practical interface that can be easily customized for the different coupled problems. Moreover, the co-simulation approach enables advantageous features. For instance, the solvers can exchange their data without having a duplication overhead, which results in saved computational time. Furthermore, as the method splits the code into functional modules (such as the data transfer, coupling operations, relaxation methods, and predictor modules), switching a given module can be done in a not intrusive way through an external configuration script. A first application of this co-simulation approach to FEM-DEM coupled analyses can be found in [53]. In this work, the co-simulation method is enhanced and extended by developing a new PFEM-FEM coupled strategy for FSI problems. As in [52,53], the coupled method is implemented in the open-source Kratos Multiphysics code.
In conclusion, this work proposes a unified open-source platform for the numerical simulation of the impact of destructive flows, such as rockslide or tsunami, against different types of retaining structures, from rock-fall cable nets to thin walls and membranes. This is purchased via an efficient co-simulation approach that couples particle methods, such as the DEM and the PFEM, with the FEM in a fully partitioned and modular way.
To the best of our knowledge, this modular approach represents the first example of a unified and open-source framework for this kind of FSI and rock/soil-structure interaction problems, which allows using such a variety of FEM structural models and particle-based methods. Indeed, depending on the problem to simulate, for the structures, we can use different constitutive models (e.g., elastic, elastoplastic, damage), solution algorithms (explicit as implicit ones), and elemental models (e.g., for solid elements, membranes, or nets), and for geophysical flow simulation, both continuum (PFEM) and discrete (DEM) particle-based methods can be employed. The modularity and generality of this uni-fied framework allow for its applicability to a large variety of relevant problems in the civil engineering field.
It is also worth remarking that in the literature there are already several efficient models that solve FSI problems that use particle methods. However, either the scope of these methods is more limited than the one proposed in this paper or, although they have a similar generality, they use commercial software and not open-source codes. For example, remarkable partitioned approaches using PFEM and FEM for FSI analysis can be found in [10,40]. [40] proposed a fully partitioned approach for FSI problems which uses an explicit PFEM for the fluid domain and an explicit FEM for the structure solved via a commercial software. The use of a commercial tool, on the one hand, allows them to model a large range of engineering problems (see, for example, the application of the method to airbag simulation in [41]), but on the other hand, it limits the widespread usability of the code. [10] proposed a fully partitioned Lagrangian framework for the solution of challenging FSI problems, using the PFEM for the fluid and the FEM for the solid. This method is based on two different solvers that are coupled through the CUPyDO framework. Differently from [10], the proposed approach takes advantage of the analogous construction of the solvers in the same framework (Kratos Multiphysics), avoiding information overhead and data duplication at the interface. In addition, we remark that the partitioned approaches proposed by [40] and [10] can only take into account the PFEM-FEM interaction and are not thought for DEM-FEM interaction, as the proposed algorithm also does.
Only very recently, partitioned approaches considering the simultaneous presence of the three methods, i.e., FEM, PFEM, and DEM, have been proposed. This is the case of the works of [11] and [48]. In particular, [11] simulated the failure of structures under the impact of tsunami waves, and [48] applied the same method to large-scale civil engineering problems. Although the present work takes some elements from [11] (as will be explained in Sect. 5.3), the methods differ for the use they make of the DEM solver. In the works of [11] and [48], the DEM is used as an auxiliary technology to model fracturing solids and mutual solid-solid interactions but not to represent a granular material avalanche, as it is done in this work.
In the literature can also be found very interesting approaches simulating the impact against retaining structures of debris flows modeled with coupled CFD-DEM, such as in [36] and [37]. These methods differ from the one presented in this work for the type of the CFD solver. For example, in [36], the CFD solution is obtained with a lattice Boltzmann method (LBM) strategy, while, in [37], following [67], a finite volume solver was used for the CFD solution. It is also important to remark that, differently from these methods, in our approach, the interaction between the DEM and the PFEM (that represents our CFD model) cannot be considered and is left for future work.
Besides giving a detailed description of the coupled solver and providing practical and useful information for new potential users, this work shows an accurate validation of the method for representative numerical tests in the field of FSI and rock-structure interaction analyses.
The paper is structured as follows: • Section 2 describes the discrete element method.
• Section 3 discusses the particle finite element method and its solution algorithm. • Section 4 describes briefly the finite element method formulation. • Section 5 delves into the co-simulation approach, making emphasis on the new parts developed in this work. • Section 6 presents five numerical examples in which the proposed approach is tested. • Section 7 presents a conclusion to this work and an outlook for future research.

The discrete element method
Unlike the PFEM, which belongs to the group of continuumbased particle methods, the DEM considers the motion and interaction of discrete particles [15]. The DEM is a particle method that decomposes the body into discrete particles. For this reason, the method is wellsuited for modeling granular materials at very different scales, from powder simulation to rockfalls, e.g. see [63] and [68], respectively.
For simplicity and computational cost considerations, generally spherical shapes are considered for the DEM particles. Nevertheless, this simplification does not always give a faithfully representation of the physical problem. In the specific case of rockfalls simulations, using spherical particles can drastically change the behavior of the solid debris and would make necessary a complex calibration of friction and contact parameters in order to have a faithful representation of the actual physical behavior. For this reason, in this work, arbitrary shapes for the DEM particles are considered. This is done using the algorithm of DEM clustering presented in [33] and discussed in detail in [6,7], where a free-to-use online tool is also provided [5].
The general workflow of a DEM-FEM simulation can be separated in: 1. Contact detection, 2. Force evaluation, 3. Integration of motion, and is elaborated on in the following subsections.

Contact detection
In contrast to using polyhedra, the usage of sphere clusters results in a reduced computational time when calculating overlap [39]. With spheres or clusters of spheres, only sphere-sphere, sphere-line, sphere-vertex, and spheresurface interactions need to be considered. We remark that these are simple operations in which only the smallest distance and the respective sphere radius are compared, which has been investigated in [50,51]. The aforementioned work additionally introduces the double hierarchy method, which describes an efficient way of handling various contact partners at the same time.
The sphere-line contact procedure is exemplary depicted in Fig. 1a. A sphere with the center C i and a corresponding radius R i is in contact with an arbitrary geometric object, in this case a line, as soon as the shortest distance is smaller than the radius d e < R i . More details about the computation of d e for different geometric entities, such as vertices, lines, and surfaces, are provided in [51].

Contact forces
After detecting the contact, the respective interaction forces must be evaluated with the chosen contact model. In the examples considered in this work, a Hertz-Mindlin springdashpot model (HM+D), described in [14], is used. The corresponding rheological model and its parameters are provided in Fig. 1b, where • k n , k t : Normal and tangential spring stiffness.
• c n , c t : Normal and tangential damping coefficients.
A detailed investigation of the HM+D model and alternative contact models can be found in [14,57,58,62], and the derivation of the contact forces F c can be found in [50,53].

Integration of motion
After the evaluation of the contact forces, the integration of motion is performed with an explicit time integration scheme. Following [39], a velocity-verlet (central differences) scheme is used to integrate the translational motion, while a quaternion-based method is employed for the integration of the rotational motion of arbitrarily shaped objects [25,30].
Using the assembled forces and torques, the translational and rotational motion is integrated according to Newton's second law. While the mass m relates the translational accelerationü to the forces F, the inertia tensor I is used to  (2)

Cluster creation
The creation of the DE clusters is depicted in Fig. 2. First, a CAD model of the desired shape is created. Subsequently, the CAD model is meshed with tetrahedral elements. The mesh is used to calculate mass and inertia. With the help of the sphere-tree algorithm presented in [6,7], the particle cluster is finally produced.

The particle finite element method
The PFEM is a Lagrangian numerical method that considers the nodes of a given mesh as particles of the fluid, which are free to move across the domain [46]. The method is well suited for free-surface flows, like dam breaks and channel flows. A comprehensive introduction of the method has been given by [12]. In the following, only the most important definitions will be recalled. The reader will be referred to the relevant sources when needed. For the following, consider a continuum domain f ⊂ R n ; n = 2 or 3 in a time interval [0, T ]. The Navier-Stokes equations govern the motion of the fluid, and defining the stress σ f = σ f (x, t) and the velocity v f = v f (x, t), (where x are the current spatial coordinates), the momentum balance and mass conservation equations read: where ρ f (x) and κ f (x) represent the fluid density and the bulk modulus, respectively, and b(x, t) the external body forces per unit mass. It must be noted that the convective term disappears because of the Lagrangian description of the fluid.
In this approach, Eqs. (3) and (4) are solved considering a slight compressibility of the fluid material [47]. Furthermore, the PFEM can consider real compressible fluids. The additional requirement is to add an equation of state and an energy conservation equation [12].
The principle of virtual work is also applied to the above equations for the computation of the solution. Following the standard PFEM, linear elements are considered and the same interpolation order is used for both the velocity and pressure degrees of freedom. To circumvent the unfulfilment of the infsup compatibility condition [12,19], the problem is stabilized using the finite increment calculus (FIC) method [47]. More details about the stabilization procedure can be found in the referenced work.

PFEM steps
A short summary of the PFEM solution steps is provided below. For a more detailed description of the method, the reader is referred to, e.g., [9] and [12]. First, consider an initial domain as depicted in Fig. 3, which consists of a fluid subdomain (which in turn may be composed by separated parts), some fixed boundaries, and a solid subdomain (shown in gray). A generic cloud of points C is also defined, which contains the set of all nodes which are part of either the solid or the fluid. Additionally, the mesh M is also defined.
The steps are as follows: 1. Fill the domain with the cloud of points n C. The n refers to the current time step t n . Note that the nodes of the fluid subdomain are also referred to as "particles". 2. Generate the finite element mesh n M. In this work, the Delaunay tessellation is used [12]. 3. Identify the external and internal boundaries for both fluid and solid domains. For this identification, the Alpha-Shape method [20] is used. This step is very important because some boundaries get severely distorted and thus need to be adequately reshaped or even removed. 4. Solve the Lagrangian form of the governing equations. 5. Update the node coordinates for the new set n+1 C, where n + 1 denotes the next time step. 6. Repeat steps 2 through 5 for the next time step.
Some remarks are necessary: • In this work, the remeshing is done at every time step. This leads to a higher mesh quality but at an increased computational cost. Nevertheless, as noted by [9], the remeshing computational time in the PFEM is not elevated if compared to the cost of the other time-step solution operations, around the 16% for a 800,000 nodes mesh according to [9]. • In step 2, the mesh elements and connectivities are redefined rather than being completely remeshed. This is because the nodes of the previous mesh are retained. For more details, the reader is referred again to [12].

The finite element method
The FEM is used to discretize the structural domains (the net structures and the flexible walls and membranes) in the present work. For this, consider an updated continuum domain s ⊂ R n ; n = 2 or 3 in a time interval [0, T ].
The momentum equation can be written as [26]: where σ is the Cauchy stress tensor, b the body forces per unit volume, ρ s the density of the body andü its acceleration. Furthermore, consider the body to be subjected to the initial conditions: where u,û 0 , v andv 0 are the displacement, its initial conditions, the velocities and its initial conditions, respectively.
The following boundary conditions must be applied to properly close the problem: In the first part, ∂ u is the Dirichlet boundary where displacements are applied. On the other part, ∂ σ is the Neumann boundary where the loads (surface tractions) are applied.
These equations are solved by means of the principle of virtual work, which applies a virtual displacement to Eq. (5).
Finally, this is solved numerically via a Newton iterative technique, which operates on a linearized version of the virtual work form of Eq. (5), as detailed in [3].

Neo-Hookean hyperelastic material model
In this work, the main constitutive laws used are the linear elastic and the Neo-Hookean hyperelastic material models, the latter being discussed in the following. This model is described by the strain-energy functional nh . A formulation provided by [3] is: where λ 0 and μ 0 are the Lamé parameters from the linear theory, C = F T F is the right Cauchy-Green tensor, J = det F and F is the deformation gradient tensor.
It is important to mention that this formulation is based on an isotropic hyperelastic material model and is valid for large deformations in a nonlinear elastic domain [3].
We note that the co-simulation approach does not impose any limitation to the use of different constitutive models. In the specific case of the computational framework of Kratos-Multiphysics, any other material law available in the constitutive model library (called ConstitutiveLawsApplication) could be used. This is one of the main advantages of the proposed co-simulation approach.

Plate in membrane action
Following by [56], membrane elements are used to model the behavior of the rock-fall cable nets under the impact of solid rocks simulated with the DEM. Moreover, the same membrane elements are also used to model flexible walls subjected to water flow impact solved with the PFEM. A detailed derivation of the membrane element and its corresponding application to model protection structures can be studied in [55]. In this work, we used the same model as in the aforementioned work, but using an isotropic elastic consistent linearized tangent modulus.
Especially for light-weight rockfall protection nets, which naturally cannot carry in-plane compression forces, special treatment of the stress configuration is needed. Following [45] the tension field (TF) theory is applied to the net models of rockfall protection structures. The TF analyzes the current principal strains and stresses and subsequently transforms the material tangent modulus, resulting in a compression-free stress state. Successful application of this theory to rockfall simulations can be found in [54,55].

Co-simulation
After presenting the individual solvers, we describe here the co-simulation approach that is used to couple the physically independent systems. A modular structure is used to combine the particle methods and the FEM in a non-intrusive way. The final result is a unified platform for simulating challenging engineering problems, such as the impact of rapid fluid flow and rockslides on retaining structures.

Modular approach
The modular approach developed in this work is based on the first implementations done in Kratos Multiphysics environment and presented in [8] and [53]. Further details of Kratos can be found in [16,17,38]. In this work, we expand these previous contributions by coupling the PFEM with the FEM in the same computational framework.
The main idea of the co-simulation approach is to couple different solvers in a non-intrusive manner and join them into a global method. In other words, the independent solvers, i.e. the particle methods and the FEM presented in the previous sections, are not modified depending on the physical problem to solve. In fact, the distinct solvers are treated as black boxes and communicate each other through an external structure: the solver wrapper. A graphical representation of the co-simulation approach is provided in Fig. 4. Several modules, i.e., coupling operations, connect the different solvers. The most important of these, e.g. the convergence loop, is described in the following sections. Details on other secondary modules available in the co-simulation approach, such as the predictor and the IO operations, are explained in more detail in [8].

Structure-particle equilibrium
For an accurate solution of the coupled problem, both particle methods (DEM and PFEM) and the FEM must be at equilibrium at each computational step. This is guaranteed by a suitable definition of the interface conditions. For the analysis of the problem, consider Fig. 5, which depicts a general particle domain p in contact with a structural domain  The data transfer between the particle methods and the FEM occurs at the particle-solid interface . There, the interaction forces arising from the particle solvers and the FEM must be equal and opposite. In equation form, this is equivalent to: where both f , p and f ,s are the forces computed at the interface and arising from the particle and structure domains, respectively.
A Dirichlet-Neumann coupling scheme is used for the partitioned strategy. This method assumes that the structure dominates the interaction [32]. Based on this equilibrium principle, the particle domain (either DEM or PFEM) is solved and the computed forces are transferred to the FEM domain. Afterwards, the FEM domain is solved and the velocities and displacements (and contact forces, for the DEM particle case) are transferred back to the particle domain to update the position of the interface.
In this work, the aforementioned data transfer takes place at every time step for the entire simulated time. As mentioned earlier, this work uses conforming meshes on both DEM to FEM and PFEM to FEM since it reduces the mapping effort to a direct transfer of the information. If this is not the case, only the mapping method needs to be changed within the setup script, such as the nearest element or barycentric mapping (see [8,65]).
Due to the distinct intrinsic nature of the two particle methods (the DEM is a discrete method, whereas the PFEM is continuous), in the DEM and the PFEM, the interaction forces are computed in different ways. In the following subsections, the algorithm of computing the interaction forces coming from the two particle methods is described separately.

DEM to FEM interaction forces
One of the main challenges of the DEM is the transmission of the coupling forces onto the FEM for the interaction. For this purpose, the Direct Interpolation Method [50] is used here, and a brief overview is given below.
Consider the schematic representation of a solid-structure interaction problem shown in Fig. 6a. The DEM particle, p , impacts a flexible net structure, s , whose boundaries are named .
For the computation of the contact forces, we use the steps described in Sects. 2.1, 2.2 and 2.3. Then, the transmission of the forces from the DEM to the FEM takes place by applying the direct interpolation method. This technique assumes that the distributed pressure (Fig. 6b) can be represented as a concentrated force at the contact point of the sphere with the interface. Then, the point forces acting through the element nodes can be computed as [50]: where N i is the shape function value at the contact point x C P per element node, and f p are the contact forces computed as discussed in Sect. 2.2. The result after performing this process is graphically shown in Fig. 6c. The aforementioned steps are summarized via pseudocode in Algorithm 1.

PFEM to FEM interaction forces
One of the main advantages of the PFEM is the easiness of its coupling with the FEM for solving FSI problems. The PFEM allows for a natural recognition of solid interfaces through the remeshing step. In the standard PFEM algorithm, the fluid and the solid elements at the interface share their boundary nodes giving rise to a conforming-mesh FSI for particle in particle_set do Loop through all sphere-particles in the particle set 3: if particle is close to then 4: Search nearest neighbors and find contact (See Sec. 2.1) 5: Compute contact_force = f , particle using the HM+D contact model (see Sec algorithm. This body-fitted method allows an easy transmission of boundary conditions at the fluid-solid interface, as it is explained below. On the other hand, having a conforming mesh FSI method can be limiting in some cases. For example, it can be problematic for thin solid interfaces embedded in a fluid or to manage very different discretization sizes of the fluid and the solid domains. In these cases, a non-conforming mesh algorithm may be preferable. Examples of these types of FSI strategies are more common in the so-called PFEM-2 ( [28]) than in the standard PFEM. An example of non-conforming mesh strategy in a PFEM-2 framework can be found in [2].
Consider the schematic representation of a FSI problem given in Fig. 7a. The PFEM domain, p , impacts a flexible structure, s , whose contours is named .
At each interface node, the impact forces are evaluated as: where n nodes is the total number of interface nodes of the element, p i is the nodal pressure, A e is the fluid element surface area (in 2D problems the element length is multiplied by an additional width parameter), N i is the element shape function value evaluated at the node and n e is the element unit normal (assumed to point towards the fluid). This is shown in Fig. 7b, c. Note that for impact problems, as those analyzed in this work, the tangential components of the shear stress tensor are usually negligible compared to the pressure normal to the face and thus can be discarded. Instead, in a general case, the normal projection of the whole Cauchy stress tensor must be considered as well [11].
An algorithm for the present case in pseudo-code is given in Algorithm 2.

Strong coupling approach
It is well known that the staggered, i.e. not iterative, solution of coupled problems may lead to a lack of accuracy and numerical instabilities of the global solution. In the specific case of DEM-FEM and PFEM-FEM coupling, the drawbacks of weakly coupled schemes have been described in [66] and [29], respectively. In this work, to avoid the mentioned numerical inconveniences, a strongly coupled partitioned scheme is used. Obtain element information: nodes, area, unit_area 4: for node in element_nodes do Loop through all the nodes in the current element 5: Get pressure_node from particle solver 6: Compute force_node = f , particle using Eq. 13 7: Transmit back the force_node values 8: end for 9: end for At each time step, the PFEM/DEM and the FEM are solved in a Gauss-Seidel iterative loop until a certain tolerance defined at the interface is reached. In particular, the convergence is considered reached when the following condition is fulfilled (see [34]): where n eq are the number of degrees of freedom of the interface, k r is the residual at the interface and is the user-defined convergence tolerance. In this work, a tolerance of = 1 × 10 −6 has been used, unless otherwise specified. As noted in [53], the interface tolerance must be larger than the convergence tolerances of the individual coupled solvers, otherwise the convergence criteria will be impossible to fulfill.
The residua of Eq. (14) are computed as the solutions difference at the structure interface between two consecutive nonlinear iterations k − 1 and k, as follows: Following [11], the relaxation is performed with respect to the residuals (16): being ω an appropriate relaxation factor ω. Furthermore, from [53], the relaxation can be performed for different variables. In the DEM-FEM case, it can be done with displacements/velocities or with contact forces. For the PFEM-FEM case, the preferred relaxation is done with the velocity only, as it is an independent variable of the PFEM solver.
In order to achieve an optimal relaxation, the Aitken method is used [31]. This method optimizes the relaxation factor ω in every iteration with respect to the current residuum s Fig. 8 Aitken relaxation diagram, adapted from [11] k r and the previous residuum k−1 r . Using the formula from [34]: The relaxation process is schematically shown in Fig. 8, for the PFEM-FEM coupling.
All the residuals are evaluated at the interface . Moreover, to keep the residual computations consistent throughout the iterations, a process to set the velocity values fixed is required for the PFEM-FEM coupling.
The reason is that the PFEM solver computes the nodal displacements by means of a Crank-Nicolson time integration scheme (using the velocity as the independent variable), whereas the FEM solver computes the displacements with an implicit Newmark-beta time integration scheme. Due to these different numerical strategies, a bias in the interface position will inevitably arise. To avoid this issue, the fixation process is implemented in the domain that does not dominate the interaction (in this case, the particle solver). The way this is done is by keeping the interface velocity values unchanged within that iteration. The displacement of the interface also remains unchanged during the solver iterations and is updated only after convergence is achieved. Further details of this fixation process can be found in Appendix A.1.

Partitioned solution algorithm
Considering all the previous discussion, an algorithm in pseudo-code is given below.
A depiction of the strong coupling algorithm is provided in Fig. 9. Algorithm 3 Staggered strong coupling 1: Initialize t n = t n−1 + t 2: for k in total_iterations do 3: Reset particle domain kinematic values ( n u p , n v p , n a p ) (See Appendix A.2) 4: Set v ,n = 0 (See Appendix A.1) 5: Solve particle domain p 6: Set v ,n = v ,n−1 7: Compute f , particle and map them to s 8: Solve FEM domain s . 9: Map n u and n v on from s to p 10: if not interface_residual ≥ tolerance_interface then Check convergence of velocities at 11: Relax the velocities at the interface via Aitken relaxation (Eq. 19) 12: Increment k and proceed from line 3 13: else 14: Update t n−1 = t n and continue to next time step 15: end if 16: end for

Results
In this section, we present five numerical tests to show the accuracy and robustness of the proposed approach. For each particle methods, i.e. the DEM and the PFEM, two tests are presented. In both cases, the first numerical test is used for validation purposes, while the second one is presented to show the applicability of the method to real cases and to provide benchmark solutions useful for future comparisons with other numerical methods. The final test is meant both to show the generality of the proposed modular approach and to provide an academic benchmark that can be useful for the validation of other solvers combining FEM with particle methods.

DEM benchmark test
This first benchmark test is used to validate the DEM-FEM coupling. As shown in Fig. 10a, a perfect sphere hits a single-span beam vertically with a given initial velocity. A coefficient of restitution of 1.0 reflects a perfectly elastic impact. To resolve the very short contact time of ≈ 0.16ms a time step of 5 10 −8 s is chosen. The structure is discretized using 60 linear beam elements.
The comparison values were taken from [50], and an analytical solution is given in [42]. The graphs of Fig. 10b show a very good agreement of the simulation data with the reference values, proving the correctness of the presented coupling DEM-FEM algorithm.

DEM large-scale civil engineering structure
This test is presented to show a practical application of the FEM-DEM coupled method to a realistic civil engineering problem. In this case, we simulate the impact of some boulders on a large-scale rockfall protection system. In cooperation with the Swiss company Geobrugg, 1 a recently published paper [55] discusses the applicability of DEM-FEM coupling to these kind of attenuator civil engineering structures. To demonstrate the generic application of the herein described coupling methodology, the aforementioned experiment is modified and subsequently simulated. While the structural setup is maintained as in [55], more rocks are simulated with varying impact positions and velocities.  The initial setup is presented in Fig. 11a, and a corresponding photograph is shown in Fig. 11b. Attenuator structures guide the impacting rocks and lead them to a safe zone. As they do not catch the rock, attenuator structures normally do not need to be replaced after each impact. They act like curtains and present a challenging FEM modeling. A variety of different structural models and material laws is necessary, as it is discussed in detail in [55].
In the present impact scenario, the rock blocks of approximately 0.75×0.51×0. 48[m 3 ] impact the attenuator structure with the given initial values.
The irregular rock shape is efficiently modeled via sphere clusters. Further information is omitted at this point, and the reader is redirected to the aforementioned publication for more information about structural parameters, modelling challenges and particle parameters. Figure 12 shows snapshots in time of the deformation pattern of the net under the impact of the solid block.
With the aim of providing useful data for future comparisons with other computational methods, we plot the time evolution of the path and velocities of the rocks in Figs. 13 and 14, respectively.

Collapse water column against elastic structure
In order to test the PFEM part of the method, the benchmark simulation of a water column hitting an elastic structure is used. This benchmark was originally proposed by [64]. The initial setup is shown in Fig. 15, and the material parameters are shown in Table 3. The fluid is discretized with 13484 triangular elements, the solid with 318, with a 2-mm average element size for both domains. The time step used is t = 0.0005s. The tolerance for the convergence check described in Sect. 5.3 is = 1 × 10 −5 .
The results of the benchmark using the proposed strong coupling algorithm for FSI are shown in Fig. 16a. They are compared against the available literature, i.e. [10,11,40,49]. Here it can be seen that the present method matches very well the literature, especially within the first half of the simulation. On the second half of it, this method also lies in between the shown results, tracking the general behavior well. As mentioned by [22], the differences between the different approaches after t = 0.5s are influenced by randomness as well: a slight change in the first part of the simulation may produce large variations in the results on the second part.
To gain more insight on the convergence acceleration, in Fig. 16b, the number of interface iterations per time step is Fig. 15 Initial configuration of the collapse water column against an elastic structure. Adapted from [11]    shown. It can be seen that the average iteration number lies between 8 and 9, spiking at the time of the first contact between the fluid and solid (at around t = 0.118s). There are a limited number of iterations in which the maximum iteration number (in this case 20) was reached, but in general the method was able to lie within this threshold. Finally, in Fig. 17a through f, several snapshots of the test are shown. There is good agreement with [21].

3D test: water slide against a flexible membrane
In order to demonstrate the capabilities of the method for realistic engineering problems, a 3D computationally demanding test was performed. The initial setup is shown in Fig. 18. It consists of a rigid container that has a flexible vertical membrane. The slanted container causes the motion of the water and its impact against the membrane, which in turn deforms and contains the fluid motion. This can be viewed as the generalization of the previous example to the 3D analysis. A further source of complexity is given by the thinness of the retaining structure whose thickness is only 1cm.
The material parameters are shown in Table 4. The fluid is discretized with 430,535 elements, the solid with 4,624, with a 2.0 cm average element size for both domains. The tolerance for the FSI iterative loop is = 1 × 10 −4 . The time step used is t = 0.001 s. Rayleigh damping [1] was used to avoid numerical noise in the simulation. For the solid, we use a St. Venant-Kirchhoff hyper-elastic material law, and for the liquid a Newtonian 2D law.
A graph showing the time evolution of the displacement of the membrane is shown in Fig. 19a. The node chosen for the analysis is marked with an A, and located +0.3 [m] above from the origin in the y-axis, as shown in Fig. 18. Also, the iterations per time step are shown in 19b.
From Fig. 19a, it can be seen that as water hits the membrane at around t = 0.420s, the membrane deforms accordingly until it stabilizes.
The largest displacement component (around 16 cm) is in X-direction, and the second largest one (around 4 cm) is in the Y-direction. This large difference was expected as the potential energy of the water body is largely converted to kinetic energy in the flow direction (X-direction) when it hits the membrane. The Z-component is close to zero, and the deviations can be explained due to randomness in the directions in which the particles move, eventually having a small imbalance and inducing some deformation in this direction.
A 3D visualization of the results is provided in Fig. 20a-e. With this example, it has been demonstrated that the proposed method can solve 3D problems as well.

2D Benchmark: water and rock long slide against a flexible barrier
In order to demonstrate the modular capabilities and ease of use of the method, a 2D setup of a representative landslide/structure interaction problem is proposed. In order to show the generality of the method, two different types of sliding materials and consequent flow regimes are considered, namely a fluid and a rock flow. The initial setup is shown in Fig. 21. The material parameters are given in Table 5. The DEM rock part is discretized with 87 elements, and a weak coupling scheme is used for the DEM-FEM examples. The PFEM fluid part is discretized with 895 elements and a Newtonian model is used. The tolerance for the strong coupling iterative loop is = 1 × 10 −4 . In both DEM-FEM and PFEM-FEM simula- (g) Velocity and displacement scale. tions, the solid structure is discretized with 40 FEM elements, for an average element size of 5 cm, Rayleigh damping is used to avoid numerical noise in the simulation, and time step duration is set as t = 0.001s. In Fig. 22a and b, we show the final position of the flexible barrier (i.e. taken when the rock and fluid flow are practically at rest), obtained with the three different solid models for the DEM-FEM and PFEM-FEM simulations, respectively. In both cases, for the elastic models, the St. Venant-Kirchhoff material law yields a smaller displacement than the Neo-Hookean law. However, while for the elastic models the rock slide gives sensibly higher displacements than the fluid flow, for the elastic-plastic models the results are similar. These different results and behaviors between the elastic and the elastoplastic models are totally reasonable and are due to the different peak impact forces given by the fluid and the rock flows. In fact, the fluid flow reaches the barrier at higher velocities than the rock flow and this translates into a transfer of a bigger momentum to the structure that induces larger plastic deformations.
This behavior is clearly shown in the representative snapshots collected in Fig. 23. The second-row pictures c, d clearly show that the fluid flow impacts at highest velocities the barrier (the rock flow reaches the obstacle after around 2.13 s at around 7 m/s, whereas the fluid flow after only 1.3 s at around 8 m/s), and this produces a larger displace-

Conclusions
We presented an advanced modular approach to couple the finite element method (FEM) with particle methods, such as the discrete element method (DEM) and the particle finite  Each computational method is used to model different physics. The FEM is employed for the solid structures: flexible membranes, nets and elastic barriers. The PFEM is used for modeling free-surface fluid flows in the presence of large topological changes. Finally, the DEM is used for rockslides simulation. To obtain a faithful representation of the actual shape of the single boulders, a clustering procedure is used in the DEM.
The distinct methods are coupled with an iterative Dirichlet-Neumann scheme. An Aitken relaxation technique is also used to improve the convergence performance of the strongly coupled algorithm.  Five numerical tests have been presented to show the accuracy of the method and its potential in the mentioned field of application. For each particle method, first, a benchmark test has been presented to prove the agreement of the proposed method to the results of the literature. Then, a more computationally demanding test has been proposed to show the applicability of the method to real-world problems. The presented numerical results have clearly shown accuracy and robustness of the method. A final test reproducing both a rock and a water flow against a flexible membrane has been proposed to demonstrate the generality of the method and to provide a benchmark problem that can be used for useful comparisons with other coupled methods. In future works, more complex constitutive models can be considered, such as non-Newtonian models in the PFEM or hyperelastic models for membranes in the FEM.   Code availability The software used for the computations is [16,17,38]. The current version is to be found as a repository in [18] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A.1 Solver wrapper fixation of velocity
As introduced in Sect. 5.3, a fixation process for the velocity in the PFEM solver is used to keep the interface displacement computation consistent throughout the iterative solution process. A deeper explanation of why this is necessary for the PFEM-FEM coupling follows. The PFEM-FEM fully Lagrangian strong coupling is based on mapping the force f (due to the pressure p) of the fluid onto the structure. As a result, the structure experiences a displacement u with a certain velocity v, which is then mapped back again to the fluid in an iterative process. As mentioned in Sect. 5.3, the PFEM uses a Crank-Nicolson time integration, whereas the FEM uses an implicit Newmark-beta time integration for the displacement calculation. Then, to avoid the convergence being affected by the difference in the interface displacement computation schemes, the PFEM displacement is fixed so that the FEM displacement dominates the interaction.
An algorithm in pseudo-code form for the PFEM case is given in Algorithm 4 to achieve the fixation.
The idea of this algorithm is to obtain and store the interface nodes. Then, all the nodes at the interface are looped through. After that, each associated equation for the velocity Algorithm 4 Fixation of velocity 1: Obtain interface_nodes and store them in memory 2: for node in interface_nodes do 3: for equation in degree_of_freedom_associated_equations do 4: Set right hand side (RHS) of equation to 0.0 5: end for 6: end for degree of freedom is looped through, and each equation's right-hand side is set to zero. What this means is that the velocities are set to zero, so when the the assembled system of equations is solved, the interface nodes have a zero velocity, which keeps their displacement unchanged and thus achieving the fixation of the values.
As a practical implementation remark, in KRATOS, a Boolean flag IsFixed is available for such purpose, and usually the BlockBuilderAndSolver classes used in the construction of the solvers have a BuildRHS method that automatically deals with all nodes flagged as "fixed".
This fixation process must not be confused with the application of boundary conditions to the domains. For the particle solver, usually non-slip boundary conditions are defined on the structure surfaces and FSI interfaces (see Fig. 24a), which can be considered as a relative velocity of the particles with respect to the surface. In contrast, the process discussed in this section fixes the global velocity of the structure moving through space (see Fig. 24b, where this global velocity is shown as v p and v s ).

A.2 Coupling operation solver reset
To reach equilibrium in Eq. 11 for the strong coupling iteration cycles, the solvers must be able to be reset. As mentioned by [34], a solver that cannot be reset is in general not suited for strong coupling computations. Therefore, it is necessary to implement a function that resets the kinematic values (displacement, velocity and acceleration) to their last converged state.
Furthermore, the algorithm must reset the particle domain values, but leave the interface values untouched. This selective reset is necessary because when the entire domain (b) Velocities of the particle and structure in the coupled case (including the interface) is reset at each iteration, the domain state never changes and thus convergence is never achieved.
On the contrary, if nothing is reset, the pressure field of the, e.g., PFEM fluid, keeps being updated and therefore an artificial pressure increase builds up, which could potentially cause a numerical blow-up of the solver. Although the blowup example refers to the PFEM fluid solver, other numerical blow-ups could appear in any kind of particle solver. The algorithm provided here is still applicable to such solvers.
Finally, Algorithm 5 is provided in pseudo-code form for this purpose.
Algorithm 5 Solver reset 1: for node in particle_domain_nodes do Modify only particle nodes 2: Obtain ( n u p ), ( n v p ) and ( n a p ) 3: Obtain ( n−1 u p ), ( n−1 v p ) and ( n−1 a p ) 4: Set node_coordinates = initial_position + n−1 u p Update nodal coordinates 5: Set n u p = n−1 u p 6: Set n v p = n−1 v p 7: Set n a p = n−1 a p 8: end for