A Partitioned Scheme for Coupling of FEM and DEM Simulations of Granular Materials

In this contribution, a simple uniaxial compaction of granular materials is considered in preparation for studying the collision behavior of double‐hull vessels filled with granular materials. The top compressing plate is deformable and modeled by the finite element method (FEM). The particles modeled using the discrete element method (DEM) are compressed by the top plate. To this end, a partitioned technique to couple the continuum approach (FEM) and particle discretization (DEM) is presented. Full agreement of the computed reaction force in DEM with the pressure applied on the top plate in FEM implies accuracy of the used coupling technique.


Introduction
There is a wide range of multiphysics engineering problems that need simultaneous application of various discretization techniques and solvers. For this, it is required to solve different fields in a coupled manner. All coupling strategies for different multi-physics problems fall into two main categories of monolithic [1,2] and partitioned [3,4] approaches. In the monolithic approach, all governing equations of different fields are solved together at the same time. On the contrary, the partitioned approach splits the entire physical problem into sub-problems. These sub-problems are solved separately, and the field quantities are exchanged at the interface [3].
In this study, an analysis of uniaxial compaction of Poraver expanded glass granules using the FEM and the DEM is presented. The granules are studied using the DEM, while the top compressing plate is considered to be deformable and is discretized with the FEM. To couple the two solvers, a partitioned solution technique is applied using the flexible and generic C++ framework comana [3]. At the end, the accuracy of this coupling technique is evaluated by a comparison of the pressure applied in the FEM solver and the reaction forces computed by the DEM solver.

Uniaxial compaction
As stated earlier, to study the collision behavior of double-hull vessels filled with expanded glass granules, a simple uniaxial compressing model as a starting point is investigated in the current work. It should be noted that the granules play the role of filler in void spaces between two hulls of the ship and transfer the contact forces of the impact to the inner hull. They can also be crushed and absorb some of the collision energy through possible marine accidents [2]. Figure 1 illustrates the twodimensional cross section of the cubic-shape compaction setup. In this setup, a bunch of granules (with initial bulk porosity of 50% in three diameters of 4, 3.125, and 2.5 mm) are sorrounded by 5 rigid walls and a moving steel deformable top plate (with Young's modulus E = 250 GPa and Poisson's ratio ν = 0.3) with the dimension of 20×20×2 mm. All the top plate faces are fixed in X and Y directions and can be moved freely in the vertical direction (Z-axis). As Fig. 1 indicates, an uniform constant pressure of 10 kPa is applied on the upper face of the top plate.

FEM/DEM coupling
The rigid walls and the Poraver granules are modeled in the DEM code MUSEN [5]. The top plate as a punch moves down to compress the particles and is assumed to be deformable. Therefore, the top plate needs to be discretized by FEM. Here, the transient analysis of the top plate is performed using the commercial FEM solver ANSYS. The interface of the deformable top plate (FEM) (hereinafter called FEM plate) and the granules (DEM) is where a coupling scheme is required (see Fig.1). To exchange the data between the DEM and FEM, the coupling interface is described by means of a triangulated mesh applying the so-called STL format. In this regard, another top plate with same dimensions and position and discretized with triangular meshes (hereinafter called STL plate) is required to be modeled in the DEM code MUSEN.
The particles apply the contact forces to the triangle cell centers (CC) in the STL plate which define the corresponding tractions on each cell (t n ). These tractions are interpolated to obtain the traction at the quadrature points (QP) of the finite elements (I QP CC (t n )). To make this procedure possible, the FEM plate surfaces need to be discretized with surface elements (SURF154 in Ansys) in addition to its volume discretization. Then the structural solver Ansys S computes the problem and 2 of 3 Section 7: Coupled problems  evaluates the displacements (d n+1 ) at the vertices of the triangulation. In order to reduce the computational effort, subcycling was used for the particle solver MUSEN P. In subcycling, N SC particle steps are performed per each structural step [6]. Afterwards, the particle solver P computes new reaction (contact) forces on the cell centers (f n+1 ) having the updated displacements on the cell vertices, and this routine is repeated for all the time steps. The coupling algorithm is depicted in Algorithm 1 in more detail. N T S in this algorithm is the number of time steps for the structural solver and the coupling scheme. Figure 2 shows the computed reaction force on the bottom face of the top plate which clearly has converged to a constant value of 4 N. This value corresponds to the summation of all forces acting on the triangular cells by the particles. Based on the static equilibrium, the total reaction force must be equal to the applied pressure (force) on the upper face of the top plate (10 kPa × 20 mm × 20 mm = 4 N). The full agreement of the computed reaction force with the applied force implies accuracy of the employed coupling technique. To evaluate the coupling solution in a more perceptual way, a softer material like rubber with E = 0.0025 GPa and ν = 0.45 was assigned to the top plate. Figure 3 shows the deformation in both STL and FEM plates after compressing the granules.