A discretization‐convergent level‐set‐discrete‐element‐method using a continuum‐based contact formulation

The level‐set‐discrete‐element‐method (LS‐DEM) was developed to overcome the shape limitation of traditional discrete element method. LS‐DEM's shape generality relies on a node‐based surface discretization of grain boundary, and it has been used to shed new light of a variety of granular mechanics applications with realistically shaped grains and structural assemblies made of unbonded building blocks. Due to the node‐based discretization of grain boundary, the original LS‐DEM is discretization‐sensitive and it suffers from divergence of the response with discretization refinement, particularly for highly compressible problems. Previous studies have identified and addressed this issue in different ways, each with its own advantages and shortcomings. Here, we propose a methodologically‐rigorous and computationally‐efficient adapted formulation which solves LS‐DEM's discretization‐sensitivity issue. It adopts the classical contact description of continuum mechanics, wherein the contact interactions are traction‐based. We demonstrate the convergence of the adapted LS‐DEM in several highly compressible cases studies, show that it is key to correctly capturing the mechanical response, and compare it to alternative formulations.

grain breakage on crushable sand and rock assemblies, 10,11 and the formation and evolution of stress transmission paths in entangled granular assemblies. 127][18] Using calibration of the normal contact stiffness k n , it has been shown to correctly capture and predict the structural response and the failure mechanism under dynamic 15 and static 16 loads.Based on a closed-form expression relating k n and the elastic modulus E, it has been used to obtain new insights on the effects of E and the friction coefficients on slab-like TIS, 17 and on their deflection capacity. 18ne of the main shortcomings of the original LS-DEM formulation, first identified and addressed in Reference 19 and later acknowledged in References 20-23, is that the response is sensitive to the surface discretization of grain boundary.Specifically, the grains become stiffer as the number of surface nodes increases, 19,22 and, as a result, the mechanical response diverges with discretization refinement.LS-DEM's discretization sensitivity is manifested mostly in highly compressible cases, that is when the forces or pressures are large enough to induce particle deformation beyond the rigid-body limit.In contrast, when the response is governed by rearrangements and rigid-body kinematics (e.g., gravity-driven particle flow), the discretization sensitivity is likely harmless because the stiffness of the grains has very little influence on the response.One can distinguish between the two regimes through the value of the compressibility number  = √ pd∕k n , 24 where p is the pressure, and d is a characteristic dimension of the particle, and k n is the normal contact stiffness. 24he first strategy proposed in the literature to address LS-DEM's discretization sensitivity is a single-node approach. 19t is based on replacing the multiple contact interactions at a contact region between two grains * with a single contact interaction in which a single representative penetration-either the deepest penetration or the average one-is used in the contact law.0][21][22] The main problems with the single-node approach, as illustrated and discussed in detail in Reference 23, are that it is possible to have the same contact force for two very different penetration states, that the contact force is prone to discontinuities (even within a single time step), and that it does not conserve elastic energy.Another problem is that, since only one penetration is used in a contact between two grains, it fails when the grains have more than one contact region, as is typical with non-convex grains.
The recently-developed volume-interacting level-set discrete-element-method (VLS-DEM) 23 has been presented as an approach that solves LS-DEM's discretization-related issues, mainly the non-conservation of elastic energy.It abandons LS-DEM's surface discretization altogether, basing instead the contact law on the volume of penetration, that is, the intersection volume of penetrating grains (i.e., level-sets).Aside from conserving elastic energy, it is likely that this approach would indeed avoid LS-DEM's discretization-sensitivity issue, even though this has not been extensively looked into in Reference 23.The two main shortcomings of the VLS-DEM approach are that it is 2-4 times slower than LS-DEM 23 (due to the computationally expansive search of the intersection volume) and that the contact law deviates from the classical continuum contact description.This deviation leads to some internal contradictions in the formulation of the volume-type contact law and it makes it necessary to rely on assumptions regarding the volumetric normal stiffness and the normal vector, assumptions whose validity has not been tested.
The objective of this study is to present a contact formulation for LS-DEM that solves the discretization-sensitivity in a straight-forward, methodologically-rigorous, and computationally-efficient manner.This formulation is based on the classical continuum contact description and its application implies minimal, yet essential, modifications to the original LS-DEM, making it not only effective, but also elegant.
Next, in Section 2, we examine LS-DEM's discretization-sensitivity on a simple two-sphere benchmark and pinpoint the reason for it in LS-DEM's original contact formulation.We then present the proposed adapted contact formulation and show that is converges in the benchmark example.In Section 3, we examine the convergence properties of the adapted formulation for several highly-compressible cases with multiple grains, and discuss their significance.

Benchmark divergence example
We illustrate the divergence of the original LS-DEM formulation in the two-sphere compression case depicted in Figure 1A, with two identical grains.Both grains are initially fixed in space such that Grain 0's rightmost point is just touching Grain 1's leftmost one.Then, Grain 1 is displaced leftwards at a constant rate along the line connecting their centers, such that its displacement  is equal to the penetration between the grains.The response parameter of interest is the contact force P.
To examine the effect of surface discretization (SD) refinement on P, we considered four progressively refined surface discretizations † , see examples in Figure 1B, and examined the corresponding P −  curves.We found the P −  curves diverge with SD refinement.We also found that SD refinement is equivalent to increasing the normal contact stiffness parameter k n , see Figure 1D, an observation we explain later in Section 2.2.The divergent behavior observed in Figure 1C underscores the need to adapt LS-DEM so that it becomes discretization-convergent, a challenge that defines the objective of this study.

Why exactly does the original LS-DEM diverge?
The root cause of LS-DEM's divergent behavior is that the inter-grain penetration stiffness, which controls the global strength and stiffness, scales linearly with the number of surface nodes.In ordinary DEM, the grains are spherical and each contact comprises a single penetration point a, see Figure 2A.In LS-DEM, due to the surface discretization, there are many contact nodes between two contacting grains, see Figure 2B.At contact node a, the normal force acting on grain i by grain j is, similarly to ordinary DEM: where k n is the normal penalty parameter with dimensions (force/penetration), and where d j,i a and nj,i a denote the penetration depth and the unit normal vector to grain j at penetration point a, respectively.
The resultant normal force F i n acting on grain i by grain j is the vector sum of nodal forces over the N contact nodes: where N is the total number of surface nodes of grain i in contact with grain j.F i,c n is obtained by integration of the normal tractions f i,c n over the contact area s. (D) The adapted LS-DEM formulation is based on the continuum description.Each node has a tributary area A a , and the resultant force F i, * n is obtained as a sum of nodal forces F i, * n,a , analogously to the continuum-based integral.The penetrations are grossly exaggerated for illustrative purposes.
By defining the average penetration d j,i a , we find that the total grain penetration stiffness, which is defined as is linear in both N and k n .The distinction between the grain penetration stiffness k grain , which refers to interaction between grains as a whole, and the normal penalty parameter k n , which refers to contact interactions at the node level, does not exist in ordinary DEM.There, N = 1 and k grain and k n are therefore one and the same.
The linearity of k grain with k n in Equation ( 3) is reflected in Figure 1D, and is common in DEM.It is useful because it allows to relate k n to the Young's modulus, either through calibration, for example, References 7 and 25, or through a closed-form expression. 16,17In contrast, the linearity of k grain with N is unique to the original LS-DEM formulation and it is problematic because, as the surface discretization is refined, N and k grain increase proportionally, leading to the divergence observed in Figure 1.

How we eliminate LS-DEM's divergence
The key to fixing the divergence of the original LS-DEM formulation is to eliminate the discretization dependence of the penetration stiffness expressed in Equation (3).To do so, we adopt the continuum-based approach illustrated in Figure 2C, wherein the contact forces are integrals of surface tractions over the contact area s, rather than sums of nodal forces.
The dimensions of the continuum-based penalty parameter k * n are those of an elastic foundation (traction/penetration), differently from k n 's (force/penetration).This approach provides the resultant normal force F i,c n and the commensurate penetration stiffness k grain with the necessary linkage to and bounding by the finite contact area s.The original formulation expressed in Equation ( 2) lacks this bounding property.
The continuum-based resultant normal force F i,c n reads: where the superscript ,c denotes the continuum-based approach, s denotes the contact area between the grains, and the elemental dF i,c n reads: where f i,c n are the normal tractions, d j,i the normal penetrations, and nj,i the surface normals.We evaluate the integral in Equation ( 4) numerically as a sum of nodal contributions at contacting nodes, as illustrated in Figure 2D.The resultant normal force F i, * n in the adapted LS-DEM therefore reads: where the superscript , * denotes the adapted formulation.The nodal normal force F i, * n,a , which is the finite analog of dF i,c n from Equation (5), reads: where A a is the tributary area of node a, see Figure 2D ‡ .Assuming that A a is taken equal to the total surface area of the grain divided by the total number of surface nodes § , the penetration stiffness k grain, * in the adapted model then reads: The "corrective" factor A a and the elastic-foundation nature of k * n are the two keys to resolving the divergence issue.Increase in N with SD refinement is off-set ("corrected") by a proportional decrease in A a , so that the sensitivity of k grain, * on the SD becomes minor, and the global response converges upon SD refinement, as will be shown in the next section.
The tangential contact formulation requires a treatment similar to that of the normal contact formulation, for similar reasons.Comparing Equations ( 1) and ( 7), the modification of the normal nodal force in the adapted formulation amounts to replacing k n with k * n ⋅ A a .The modification of the nodal tangential force in the adapted formulation is analogous to this replacement.In the original LS-DEM (see eq. ( 15) in Reference 7), the elastic tangential force at node a, which is the basis for the commonly-used Coulomb friction law, is based on the product k s ⋅ ΔS a , where k s (force/displacement) is the tangential stiffness analogous to k n and where ΔS a is the (incremental) tangential slip vector analogous to the penetration vector d j,i a ⋅ nj,i a .Hence, the modification of the nodal shear force amounts to replacing k s with k * s ⋅ A a , where the distributed tangential stiffness k * s has dimensions (traction/penetration).The major advantage of the adapted contact formulation compared to the volume-based approach 23 and the single-node approach 19 is that it is derived directly from the well-established continuum contact description based on surface tractions.As such, it account for the spatial variation of contact tractions across the contact area, preserves the classical significance and dimensions of the penetration stiffness (traction/penetration), and obtains the normal and tangential vectors directly from the geometrical description, without any associated assumptions.All these elements make the proposed approach not only methodologically sound, but also numerically smooth, stable, and economical.

Adapted LS-DEM formulation
With the modifications to the contact formulation described above, the adapted LS-DEM's formulation is neatly obtained from the original one by replacing k n and k s in Equations ( 5) and ( 9) from Reference 7 with k * n ⋅ A a and k * s ⋅ A a , respectively.Regarding the computational aspect, it is pointed out that the additional computational cost in the adapted formulation is effectively zero, consisting merely of calculating A a once for each node in the pre-processing stage and appending it to k * n and k * s whenever the node is in contact.This speaks to the computational efficiency of our approach.

RESULTS AND DISCUSSION
The convergence properties of the adapted LS-DEM in different deformation-governed cases are demonstrated and discussed next.

Two-sphere central-compression
We begin exploring the convergence of the adapted LS-DEM formulation in the illustrative two-sphere benchmark for which the original LS-DEM diverged, see Figure 1.The force-displacement with progressively refined SD's in Figure 3A show that the adapted LS-DEM converges upon SD refinement.Further, for a given surface discretization, the response scales linearly with k * n , see Figure 3, reflecting the linearity of k grain, * with k * n in accordance with Equation (8).These results indicate that the cause of divergence in the original LS-DEM formulation has been correctly identified and eliminated as described in Section 2. (A)

F I G U E 4
Comparison between the original and the adapted LS-DEM formulation for a case of a 3 × 3 × 3 ordered-sphere assembly.

Spherical grains compression
Next, we consider an ordered assembly of 3 × 3 × 3 = 27 spherical grains with the objective of examining the convergence properties of the adapted LS-DEM formulation for assemblies with multiple grains.The examined 3 × 3 × 3 assembly is supported by fixed vertical walls and a horizontal bottom floor and subject to unidirectional compression applied by a moving upper horizontal wall, see Figure 4A.
Comparing Figure 4B with Figure 4D, we again see that the adapted LS-DEM converges with SD refinement, whereas the original LS-DEM diverges.Also, both the original and the adapted LS-DEM scale similarly (and linearly) with k n and k * n , respectively, see Figure 4C,E.These results show that the convergent behavior of the adapted formulation is not limited to the previously shown two-grain benchmark, but that it also holds for ordered assemblies of multiple spherical grains.

A slab-like assembly of interlocked blocks
The real strength of LS-DEM is its ability to deal with irregularly-shaped grains.Therefore, as a final check on the convergence properties of the adapted LS-DEM formulation, we consider a structural assembly made of such irregularly shaped grains, or blocks.Figure 5A depicts a truncated polyhedral grain, or block, of the assembly we consider.A five block interlocking formation and the entire slab-like assembly are shown in Figure 5B,C.The examined assembly was chosen for two reasons.First, its blocks have conforming planar faces with sharp edges and corners, posing an interesting challenge to the adapted LS-DEM formulation, and providing a different context to examines its convergence properties.Second, it has been studied experimentally in Reference 26, providing us with important information about its mechanical response which we will discussed later in Section 3.5.
The assembly was loaded by a spherical indenter which was displaced downwards at a constant speed, pressing against the central block with force P. Additional details regarding the configuration and the simulations can be found in Reference 16.
The P −  curves in Figure 5E-H reveal the same trends observed in the previous examples with spherical grains, namely that the adapted LS-DEM converges with SD refinement whereas the original one does not, and that both scale linearly with k n and k * n , in accordance with Equations ( 3) and ( 8).This shows that the adapted LS-DEM is not only convergent in the case of spherical grains with non-conforming contacts, but also for assemblies made of planar-face grains interacting across conforming surfaces, sharp edges and corners.

The relation between surface discretization and the level-set volume discretization
In addition to the surface discretization of LS-DEM we have focused on thus far, a uniform Cartesian volumetric grid is used in LS-DEM to discretize the level-set representation of the grains.Differently from the surface discretization, this volume discretization is not associated with convergence issues because it does not affect the total grain penetration stiffness in Equation ( 3).This is why we have not yet discussed it.The volume discretization (VD) is, however, important to discuss in the present context because its refinement, defined by the step size of the volume grid on which the level-set is evaluated, determines the accuracy with which the magnitude and the direction of the penetration at a contact node in Equations ( 1) and ( 7) are evaluated.
In the context of the original LS-DEM, it is not clear how to determine the relation between the SD refinement and the VD refinement in order to obtain consistent results.The convergence of the adapted LS-DEM established in the previous sections provides a clear-cut criterion for this determination, namely that the SD and VD should be refined such that the results do not change with further refinement of either or both together.
In this section, we illustrate a two-way SD-VD convergence study in the context of the slab-like case in Figure 5.It provides an estimate to how refined the SD should be relative to the VD in general.Differently from the previous examples, where the number of surface nodes Nnodes was taken as the SD parameter, here we choose the average distance between surface nodes as the surface discretization parameter SDP.Similarly, we choose the step size of the level-set volume grid as the volume discretization parameter VDP.The reason is that length parameters are a more meaningful measure for quantifying the relation between the SD and VD refinement than the number of nodes, which scales differently for SD and VD with the overall lineal dimensions of the grain.
Considering first the SD convergence with a constant VDP = 0.025 mm in Figure 6A, we find that SDP = 0.06 mm (with Nnodes = 68,314) is a converged SDP.Note that the larger range of SDP's (and commensurate Nnodes) in Figure 6A compared with Figure 5E allows to see the full convergence picture missing in Figure 5E, where convergence is already reached with Nnodes ≈ 17,000 (SDP ≈ 0.14 mm).Considering next the VD convergence with the converged SDP = 0.06 mm, we find that VDP = 0.025 mm is a converged VDP.Based on the converged values, we find that 0.025∕0.06= 0.416 is a suitable VDP-to-SDP ratio for convergence.This ratio provides a first-order approximation to facilitate future two-way SD-VD convergence studies.
We note that, as Nnodes grows by more than two orders of magnitude (from 420 with SDP = 0.68 mm to 68,314 with SDP = 0.06 mm), the maximal P varies by a factor of only about 2. This speaks to the small sensitivity of the adapted formulation on the SD refinement.This small sensitivity means that reasonable estimates of the converged behavior can be obtained with computationally-cheaper coarse SD's.

The physical significance of convergence
How important is it to use a converged discretization in terms of correctly capturing the mechanical response?We address this question next in the context of the slab-like assembly in Figure 5 by comparing the failure mechanism as obtained with the adapted LS-DEM to the experimentally observed one in Reference 26.
In Reference 26, the failure mechanism was described as follows: "The architectured panels failed by a progressive 'push-out' of the center block by the loading pin, while the rest of the panel remained largely intact."The failure mechanisms as obtained with the adapted LS-DEM are depicted in Figure 6C,D, through six snapshots corresponding to  = 0, 1, … , 5 mm.We find that, with the coarse discretization (SDP = 0.68 mm and VDP = 0.43 mm), the failure mechanism is qualitatively different from the experimentally observed one.The blocks stick to one another throughout most of the response and several of the central blocks end up falling-off, see Figure 6C.In contrast, Figure 6D shows that, with the converged discretization, LS-DEM captures the failure mechanism observed in the experiments as described in Reference 26.This speaks to the importance of having a converged LS-DEM formulation and working with converged discretizations for capturing the failure mechanism.

SUMMARY AND OUTLOOK
We have presented an adapted contact formulation of the LS-DEM aimed at solving the discretization sensitivity of the original formulation.Differently from existing alternative approaches, our approach is based on the well-established traction-based continuum description of contact.As such, it is mechanically sound and straightforward to understand and implement.In addition, it is numerically smooth and stable, and its computational cost is effectively zero.We have shown in a variety of highly-compressible applications (where the original LS-DEM's discretization-dependence is most problematic) that the adapted LS-DEM is fully discretization convergent.We have further shown that this convergence is key to correctly capturing the mechanical response.The proposed convergent formulation elegantly solves one of LS-DEM's major original drawbacks, paving the way for future methodological improvements.

1
Benchmark example for divergence in the original LS-DEM with surface discretization refinement under displacement-controlled loading.(A) Configuration of the two-sphere central-compression benchmark (B) two different triangle-mesh-based surface discretizations (C) force-displacement (P − ) curves with varying number of surface nodes (Nnodes) and a fixed k n = 1 GN/mm (D) P −  curves with varying k n 's and a fixed Nnodes ≈ 8000.
Illustration of differences between the contact descriptions in DEM, original LS-DEM, and adapted LS-DEM.(A) DEM contact approach: Penetration of grains i and j is defined by the position of the single penetration point a. (B) The original LS-DEM contact approach: Surface discretization nodes are indicated by black dots.Penetration d j,i a of grain j into grain i at discretization node a is shown as a dashed line.The resultant normal force F i n is obtained as the sum of nodal forces F i n,a .They are indicated by full lines and a dashed one, respectively.(C) In a continuum-based contact description, the resultant normal force

F I G U R E 3
Force-displacement curves for the two-sphere central-compression benchmark case with the adapted LS-DEM formulation: (A) P −  curves for a fixed k * n = 1 GPa/mm and five progressively refined surface discretizations.(B) P −  curves for a fixed Nnodes ≈ 8000 and varying k * n .

(
A) The assembly's configuration and boundary conditions.(B) Load-displacement P −  curves obtained with the original LS-DEM formulation with a constant k n = 1 GN/mm and varying surface discretizations.(C) P −  curves obtained with the original LS-DEM formulation with a constant Nnodes = 4000 and varying k n 's.(D) P −  curves obtained with the adapted LS-DEM formulation a constant k * n = 1 GPa/mm and varying surface discretizations.(E) P −  curves obtained with the original LS-DEM formulation with a constant Nnodes = 4000 and varying k * n 's.Each thick-lined curve represent the average of three realizations with slightly varied initial orientation of the spheres.

5
of the adapted LS-DEM in a case with planar-faced grains with conforming contacts.(A) Geometry of a single interlocked grains.(B) Interlocking pattern of five neighboring grains.(C) Entire slab-like interlocked assembly with its boundary and loading conditions.(D) Illustration of the triangle mesh surface discretization of a grain with Nnodes ≈ 17,000.(E-H) Force-displacement P −  curves of the original and adapted LS-DEM.

6
The adapted LS-DEM is fully discretization convergent.P −  curves convergence with (A) surface discretization refinement, and (B) volume discretization refinement (The VDP = 0.43 mm case has an SDP = 0.68 mm).Failure mechanism convergence from (C) coarse to (D) refined discretizations.