A medial axis based method for irregular grain shape representation in DEM simulations 

Get Complete Project Material File(s) Now! »

Computational considerations

DEM computations are performed using an explicit integration scheme, which is inherently only conditionally stable, providing that the simulation timestep Δt is sufficiently small. The critical timestepΔ tCR for continuum equations is based on the sonic speed: the elastic wave must not propagate farther than the minimum distance of integration points during one timestep. For a system of spherical discrete elements, the critical timestep can be expressed as [Zhao, 2017]:
where mi and Ki are the mass and sti↵ness of discrete element i.
The numerical cost of a simulation is inversely proportional to the timestep and since the critical timestep directly depends on the particle masses and contact sti↵nesses, the choice of these parameter values has a crucial e↵ect on the computational time.
• Contact sti↵ness: In DEM simulations, particle overlap is used as a penalty parameter to evaluate contact force, based on the contact sti↵-ness (Eq. 3.26). However, the validity of the simulations is conditioned by the rigid grain limit [Cundall and Strack, 1979; Da Cruz et al., 2005], i.e. sufficiently large contact sti↵ness to insure these overlaps are neg-ligibly small with respect to particle sizes. In YADE, contact sti↵ness is determined by the value of the radii and Young’s moduli of the two particles in contact (Eq. 3.24). Increasing Young’s modulus thus in-creases the computational accuracy, but reduces the critical timestep (Eq. 3.73) and a balance is sought between the computational accuracy and numerical cost.
• Particle masses: Particle masses are computed from the material density and particle volumes. Although material density should corre-spond to the properties of the simulated material, it can be increased in order to increase the critical timestep (Eq. 3.73) and decrease the com-putational cost. The approach is valid, providing the system remains in a quasi-static regime during the entire duration of the simulation [Sheng et al., 2004; Hagenmuller et al., 2015].

A medial axis based method for irregu-lar grain shape representation in DEM simulations

Abstract This paper describes a novel method for representing arbitrary grain shapes in Discrete Element Method (DEM) simulations. The method takes advantage of the efficient sphere contact treatment in DEM and ap-proximates the overall grain shape by combining a number of overlapping spheres. The method is based on the medial axis transformation, which de-fines the set of spheres needed for total grain reconstruction. This number of spheres is then further diminished by selecting only a subset of reconstruct-ing spheres and opting for a grain approximation rather than a full grain reconstruction. The e↵ects of the grain approximating parameters on the key geometrical features of the grains and the overall mechanical response of the granular medium are monitored by an extensive sensitivity analysis. The results of DEM quasi-static oedometric compression on a granular sample of approximated grains exhibit a high level of accuracy even for a small number of spheres.
The microstructure of a granular material, i.e. the relative arrangement and orientation of these grains, can have a crucial e↵ect on its mechanical response [Kochmanov´ and Tanaka, 2010; Zhao et al., 2016]. In turn, mi-crostructure can be strongly a↵ected by the shape of grains [Szarf et al., 2009; Albaba et al., 2015; Markou, 2016]. Nowadays the 3D shape of grains and microstructure can be directly accessed with techniques such as X-ray tomog-raphy [And`o et al., 2012; Hagenmuller et al., 2013; Schleef et al., 2014]. In this study, we will consider the case of a specific material, namely snow, for which we expect a particularly strong e↵ect of microstructure on the macroscopic behaviour [Schweizer and Jamieson, 2001; Schweizer et al., 2003]. Snow con-sists of a highly porous continous ice matrix (void ratio up to 9) that, in fast deformation regime, can be regarded as a cohesive granular material. In this regime, which is typical of snow avalanche initiation [Gaume et al., 2015, 2017b], macroscopic deformations are primarily caused by rearrangements of grains [Hagenmuller et al., 2014b] that can come in a large variety of shapes, from rounded to faceted or elongated [Fierz, 2009]. The objective of this work is to devise an efficient method to investigate the mechanical behaviour of this material through DEM simulations.
In Discrete Element Method (DEM) [Cundall and Strack, 1979; Radjai and Dubois, 2011], granular materials are usually modelled as collections of spherical particles. The deviation of the particle shape from a sphere can be taken into account indirectly through parameters such as rolling friction that simulate grain interlocking [Ai et al., 2011]. However, a priori estimating the rolling friction parameter that could capture the e↵ect of the grain shapes on the mechanical response of the granular assembly remains problematic [Wensrich and Katterfeld, 2012]. Alternatively, the exact 3D shape of grains can also be explicitly modelled in the DEM. However, assigning the shape of the grains to the discrete elements (DEs) themselves, by modelling them as polyhedron [Hart et al., 1988; Hogue, 1998; Lee et al., 2009], turns out to be very computationally expensive due to complex contact detection and contact forces calculation. A more efficient way of modelling an arbitrary grain shape in DEM is by utilising a multitude of spheres in order to capture the geometry of the grain [Kruggel-Emden et al., 2008]. A grain is thus represented by a clump of spheres that behaves as a single rigid body. This way the efficiency of the DEM in handling contacts between spheres can be preserved. The downside of the method is, since every grain is represented by a multitude of spheres, a large number of elements for contact detection, which renders the DEM simulations slow. This makes the efficiency of the grain approximating method, in terms of the number of utilized spheres, of primary importance.
Regarding the sphere-based representations of grains, a straight-forward approach was adopted by Hagenmuller et al. [2015], who modelled grains by arranging a set of smaller spherical discrete elements on an orthogonal grid along the grain boundary. This approach produces a very fine grain shape dis-cretization, but results in a very large number of discrete elements, rendering simulations very slow. Multiple authors have instead chosen the overlapping sphere approach, where the grain is represented by a set of larger spheres allowed to overlap that fill the volume of the grain. Most of the studies that adopted this approach [Price et al., 2007; Ferellec and McDowell, 2010; Amberger et al., 2012; Gao et al., 2012] are based on first populating the grain shape with a large number of spheres, and then applying a secondary algorithm that discards certain spheres in accordance with a predefined cri-terion in order to achieve a grain approximation at a moderate number of spheres. The original population of spheres in these methods can however involve redundant information [Coeurjolly et al., 2008]. In addition, some of these methods [Price et al., 2007; Ferellec and McDowell, 2010] involve a certain level of randomness in the initiation of the process. In this study, in search of an optimal grain mapping with a reduced number of spheres, we propose a method for reducing the redundancy while objectively calculating sphere positions.
Similarly to the above presented sphere overlapping methods, the ap-proach developed here is based on first creating an initial set of sphere can-didates, and then filtering this set in order to obtain a final grain approxi-mation. However, the redundancy of the initial set of spheres is significantly reduced by applying the medial axis transformation before applying the fil-tering. Thus, the initial set of spheres is diminished while still including all the spheres needed for a complete grain approximation. Although medial axis is a commonly used tool for shape approximation in computer graphics and computer vision [Bradshaw and O’Sullivan, 2004; Stolpner et al., 2012], to our knowledge it has not yet been applied in the field of granular me-chanics. In the second step, this set of spheres is then filtered in order to obtain a grain approximation at a moderate number of spheres. Compared to [Bradshaw and O’Sullivan, 2004; Stolpner et al., 2012], we employ here a simple and robust filtering criterion that allows a variety of di↵erent grain approximations depending on the choice of two approximating parameters.
The paper is organized as follows: The developed method for grain ap-proximation is presented in section 4.2 and the tools for geometrical evalua-tion of approximated grain shapes are described in section 4.3. The function-ing of the method is demonstrated in section 4.4, by applying it to a sample of snow. X-ray-tomography-derived shapes of snow grains are approximated with the developed method and the key geometrical features of the grains are studied with respect to the values of approximating parameters. Finally these results are put into perspective by simulating and studying the me-chanical response of the snow sample in quasi-static oedometric compression with respect to these grain approximating parameters.

A medial axis based grain approximating method

The input data utilized by the method is a discrete binary image of the grain assembly (in our particular case, the grains are derived by segmentation from tomographical 3D images of snow microstructure [Hagenmuller et al., 2013]). The first stage of the grain approximating process consists of performing a medial axis transformation on the grain image. This results in a set of points that represent the centres of a set of spheres necessary for a complete grain reconstruction.

The medial axis transformation

The medial axis (M A) [Blum, 1976] is a classical tool for shape recognition. It consists of all the points having more than one nearest point on the object’s boundary, thus forming the topological skeleton of the object. By combining the M A with the distance transformation (DT ), which bears the information on the Euclidian distance of each object point to the closest surface point, the information on object shape is highly reduced, yet sufficient for a complete object shape reconstruction. The object is reconstructed with a minimal necessary number of spheres by placing spheres of appropriate radii on each point of the M A (Fig. 4.1).

Table of contents :

1 Introduction 
1.1 Context
1.2 Snow avalanche formation
1.3 Snow types
1.4 The scope of the thesis
2 State of the art: snow under shear and normal loading 
2.1 Experimental investigations
2.1.1 Constitutive mechanical behaviour of snow
2.1.2 Slab-weak layer system response
2.2 Numerical modelling approaches
2.2.1 Constitutive mechanical behaviour of snow
2.2.2 Slab-weak layer system response
2.3 Scientific questions and approach
2.4 Structure of the thesis
3 Discrete element method 
3.1 YADE computing procedure
3.1.1 Collision detection
3.1.2 Strain evaluation
3.1.3 Contact force evaluation
3.1.4 Motion integration
3.2 Computational considerations
4 A medial axis based method for irregular grain shape representation in DEM simulations 
4.1 Introduction
4.2 A medial axis based grain approximating method
4.2.1 The medial axis transformation
4.2.2 The grain approximating process
4.3 Geometrical evaluation of grain approximation
4.4 Application of the developed grain approximating algorithm to a granular material
4.4.1 Geometrical sensitivity analysis
4.4.2 The e↵ect of grain approximation on the mechanical behaviour of the granular material
4.5 Conclusions
4.A Modelling grain contacts
5 Mixed-mode loading simulation of snow 
5.1 Samples
5.2 Boundary conditions and loading
5.3 Material parameters
5.3.1 Density
5.3.2 Young’s modulus
5.4 Timestep
5.5 Grain approximation
5.6 Conclusions
6 Snow failure modes under mixed loading 
6.1 Introduction
6.2 Methods
6.3 Results
6.4 Discussion and conclusions
7 A microstructural investigation of snow failure under mixedmode loading
7.1 Introduction
7.2 Discrete element simulations
7.3 Macroscopic response
7.4 Microscopic investigation
7.4.1 Contact force distribution
7.4.2 Force chain analysis
7.4.3 Force distribution within force chains
7.4.4 Force distribution around force chains
7.5 Conclusions
7.A Analysis of the spatial distribution of damage in the snow sample
7.A.1 Spatial distribution of damage at failure
7.A.2 Damage clustering
8 Conclusion and perspectives
8.1 Conclusion
8.1.1 Numerical tools for snow mechanics
8.1.2 Insights into mixed-mode loading response of snow
8.2 Perspectives

GET THE COMPLETE PROJECT