Surface reconstruction with spherical harmonics and its application for single particle crushing simulations

Deheng Wei ,Budi Zho ,Yixing Gn ,*

aSchool of Civil Engineering,The University of Sydney,Sydney,NSW,2006,Australia

b School of Civil Engineering,University College Dublin,Dublin,Ireland

Keywords:Numerical modelling Particle morphology Particle crushing/crushability Particle-scale behaviour Sands

ABSTRACT Particle morphology has great in fluence on mechanical behaviour and hydro/thermal/electrical conductivities of granular materials.Surface reconstruction and mesh generation are critical to consider realistic particle shapes in various computational simulations.This study adopts the combined finitediscrete element method(FDEM)to investigate single particle crushing behaviour.Particle shapes were reconstructed with spherical harmonic(SH)in both spherical and Cartesian coordinate systems.Furthermore,the reconstructed surface mesh qualities in two coordinate systems are investigated and compared.Although the ef ficiency of the two SH systems in reconstructing star-like shapes is nearly identical,SH in Cartesian coordinate system can reconstruct non-star-like shapes with the help of surface parameterisation.Meanwhile,a higher triangular mesh quality is generated with spherical coordinate.In single particle crushing tests,the low mesh quality produces more fluctuations on load-displacement curves.The particles with more sur ficial mesh elements tend to have a lower contact stiffness due to more contact stress concentrations induced by complexity of morphology features and more volumetric tetrahedral elements.The fracture patterns are also in fluenced by mesh quality and density,e.g.a particle with fewer mesh elements has a simpler fragmentation pattern.This study serves as an essential step towards modelling particle breakage using FDEM with surface mesh directly from SH reconstruction.?2022 Institute of Rock and Soil Mechanics,Chinese Academy of Sciences.Production and hosting by Elsevier B.V.This is an open access article under the CC BY-NC-ND license(http://creativecommons.org/licenses/by-nc-nd/4.0/).

1.Introduction

Particle morphology is one of the inherent soil characteristics determining properties of granular soils,e.g.mechanical responses(Cho et al.,2006;Yang and Luo,2015),hydro-mechanical behaviours(Gordon et al.,2004;Suh et al.,2017),and electrical/thermal conductivities(Choi et al.,2001;Friedman and Robinson,2002;Lee et al.,2017).Numerical simulations promote understanding of behaviours of granular soils using spherical packings,such as compressibility(Minh and Cheng,2013),shear strength(Cui and O’Sullivan,2006),wave propagation(Mouraille and Luding,2008)and hydraulic conductivity(Gan et al.,2013).However,numerically quantifying the relation between sand particle breakage and particle morphology still remains a challenge.

Discrete element method(DEM)is widely implemented in geotechnical engineering to simulate grain breakage with two approaches,i.e.bond method(McDowell and Harireche,2002;Cheng et al.,2003)and replacing method(de Bono and McDowell,2018;Ciantia et al.,2019).Replacing method does not satisfy mass conservation,while bond method needs careful calibration of inter-particle contact properties to simulate elastic properties,e.g.Young’s modulus and Poisson’s ratio.The combined finite-discrete element method(FDEM)is an alternative approach that can simulate continuum behaviour,fracture initiation/propagation(Munjiza,2004;Ma et al.,2016),and the large number of contacts within granular assembly.Due to the limitations of classical DEM in investigating particle crushing behaviour,FDEM is an alternative for better incorporating particle shape with the help of triangular sur ficial meshes.

Recently,micro X-ray computed tomography(CT)can obtain three-dimensional(3D)particle morphology of natural sands.The voxelised CT images with stair-step artefacts are either smoothed by marching cubes algorithms(Alshibli et al.,2015;Kong and Fonseca,2018)or mathematically reconstructed with spherical harmonics(SH)(Wei et al.,2018a).Compared with marching cubes algorithm,SH method is powerful not only in shape quanti fication and surface discretisation(Garboczi,2002),but also in generating virtual particles with speci fied shape parameters(Zhao et al.,2017;Wei et al.,2018b).Previously,the two-dimensional(2D)equivalence of SH in spherical coordinate system,Fourier analysis,is widely used in geotechnical engineering(Bowman et al.,2001;Zhou et al.,2015),while SH in Cartesian coordinate system is applied in medical images(Shen and Makedon,2006;Yotter et al.,2011).However,there are few literatures to clarify applicability and limitations of the two SH methods.

SH construction can be performed in two forms using Cartesian and spherical coordination systems,respectively.The major advantage of SH-Cartesian is the ability to reconstruct non-starshaped particles.For example,by quantifying fractured grains from single sand crushing in Zhao et al.(2015),it was found that fragments are mostly non-star-shaped.Since then,many studies started to use Cartesian-SH to reconstruct CT-scanned shapes of sands or their fragments instead of spherical-SH(Zhou et al.,2015;Su and Yan,2018;Nie et al.,2020).However,Cartesian-SH is more complicated and has more extra steps(i.e.surface parameterisation)to be conducted and more SH coef ficients.Meanwhile,when generating random shaped grains by inverse SH analysis,how to merge three sets of SH coef ficients into one value of fractal dimension needs more considerations.In addition,some key questions remain to be answered:(1)What will the non-star-shaped grain be if it is reconstructed by spherical-SH?and(2)How is the ef ficiency of two kinds of SH analysis in approximating star-shaped grains,and furthermore,mechanical-related computational modelling?To answer these questions,this study is dedicated.

This study focuses on in fluences of SH-generated mesh qualities on FDEM-simulated particle breakage behaviour.Firstly,the hybrid FDEMfor simulation of single particle crushing is simply introduced.Then the process of generating surface mesh via SH is described.Furthermore,for comprehensive comparison between two SHgenerated meshes,we reconstruct two kinds of particles with different irregularities,i.e.Leighton Buzzard sand(LBS)and their fragments(Fig.1).Subsequently,shape parameters of reconstructed surfaces in both coordinate systems are compared.Finally,FDEM simulation results demonstrate the in fluence of mesh quality on single particle crushing and are compared with those of relevant extended finite element method(XFEM)studies in Druckrey and Alshibli(2016).

2.Surface mesh reconstruction with SH

2.1.Spherical harmonics(SH)

Any spherical scalar function can be decomposed as the sum of SH:

SH reconstructs a particle surface by fitting the coordinates of discretised points on this surface.These points can be represented either in a spherical system asri(θ,φ)or in Cartesian system asvi(xi(θ,φ),yi(θ,φ),zi(θ,φ)),whereiindicates theith point.

For spherical-SH,to reconstruct a surface with SH in the spherical coordinate system is to find a centroid as a reference to calculate its distance to surface.The centroid is selected as averaged coordinates of all vertices(Zhou and Wang,2015).Once the centroid is determined,surface points must be uniquely represented by(θ,φ)which means only star-like shapes can be regenerated.Star-like shapes have at least one internal point that connects all surface point within the shape itself.If one shape does not have this property,it is non-star-like.The second step is to determine the SH coef ficients up tonmaxdegree with a cloud of surface points.A surface withNpoints results inNlinear equations with(nmax+1)2unknown coef ficients,usuallyN＞(nmax+1)2.The least-square fitting is used to find the optimised coef ficients.Whenn=0,a perfect sphere is generated,as shown in Fig.2.With the increase ofn,more and more details can be depicted.

For Cartesian-SH,surface reconstruction requires a one-to-one mapping between Cartesian and spherical coordinates,v(x(θ,φ),y(θ,φ),z(θ,φ)),called surface parameterisation.Parameterisation projects particle surface onto a unit sphere surface through a constrained optimisation process.The generated bijective mapping,between each Cartesian surface point and a pair of spherical coordinates,makes reconstructing non-star-like shapes possible(Brechbühler et al.,1995).The second step is similar to that of spherical-SH.Three systems of linear equations are established forxi(θ,φ),yi(θ,φ),andzi(θ,φ):

Fig.1.Surface meshes of an LBS particle(a)and an LBS fragment(b)with 1500 triangular elements.Hot colour indicates high curvature value,while the blue colour indicates negative curvature values.

Fig.2.Surface reconstruction at varying maximum SH degree,n max,with Cartesian-and spherical-SH for an LBS particle(a)and an LBS fragment(b).The continuum and smooth reconstructed surfaces are discretised with 320,1280 and 5120 triangular elements.Hot colour indicates high curvature value,while the blue colour indicates negative curvature values.

Three sets of coef ficients are determined,resulting in 3(nmax+1)2coef ficients.Cartesian-SH avoids the illness in reconstructing non-star-like shapes,and manually selecting the centroid in spherical-SH.Notably,Cartesian-SH does not need pre-de fined centre,which can be explicitly expressed by Eq.(2)whenn=0,which is different from spherical-SH.As demonstrated in Fig.2,a perfect ellipsoid is depicted whenn=1.

2.2.Reconstruction for sand particles and fragments

We use 3D morphologies of LBS particles and their fragments obtained from micro X-ray CT(Zhao et al.,2015;Zhao and Wang,2016).LBS particles with size equal to 1-2 mm are smooth and well-rounded resulting from the abrasion during water transportation in rivers(Krumbein,1941),while LBS fragments(0.2-1 mm)are elongated and angular(Fig.1).Similar to Ottawa sand,LBS is also a quartz sand,and its main component is silicon dioxide.LBS particles are typically star-like,while LBS fragments are likely to be non-star-like.LBS fragments are included to test the ef ficiency of two kinds of SHs in depicting non-star-like shapes.

Image processing to deal with CT images is not the focus of this study.Without loss of generality,only brief details are summarised here,and more information can be referred to our previous studies(Zhao et al.,2015;Zhao and Wang,2016).The 3D CT data can be visualised as a stack of 2D images,each with a cross-section of the scanned sample with the thickness of one voxel of high resolution of 3.3μm.The main steps include the use of a 2D median filter to reduce the noise on every slice,followed by binarisation of the filtered image for detection of the LBS particle boundaries from distinct material phases.Then,marching cubes algorithm is used to reconstruct voxelised grain shapes with 1500 triangular sur ficial elements,as shown in Fig.1.Finally,surfaces of the labelled CT images of grains or fragments are used for SH analysis.

Various sets of(θ,φ)from icosahedron-based geodesic structures are directly imported into Eqs.(1)and(2)to reconstruct irregular shapes.By repeatedly splitting triangular faces of a unit icosahedron into four similar triangles at middle points of its edges,the resulting approximated spheres have 20×4Nfaces and 10×4N+2 vertices.The 752 vertices on each particle surface are reconstructed using both Cartesian-and spherical-SH with varyingnmax.The LBS particle and fragment are reconstructed withnmax=1,5,10 and 15,as illustrated in Fig.2.We assume that the surface with 5120 triangular elements is precise enough to capture all surface features withnmaxup to 15.For LBS particles,a highernmaxleads to gradually better reconstruction of surface features with both Cartesian-and spherical-SH.The surfaces reconstructed with spherical-SH have some bumps even with the highest reconstruction degree.The more irregular fragment can result in high protrusions on the surfaces reconstructed with spherical-SH due to the local angular features violating the star-like shape rule.As shown in Fig.2,the Cartesian-SH provides a better reconstruction according to the visual inspection.To quantify shape indices of reconstructed shapes,we use four shape indices(Fig.3),i.e.aspect ratio,Ar,to compare three principal dimensions;roundness,R,to quantify the local mean curvature distribution at corners;sphericity,S,to describe how the particle is close to a sphere;and convexity,CX,to quantify the relative volume of concave features.At lownmax(i.e.nmax＜5),Cartesian-SH approximates original shapes well with relative error less than 10%except for roundness.At mediannmax(i.e.5≤nmax≤10),Cartesian-SH still works better,which can be inferred by roundness values.However,at highnmax(i.e.nmax＞10),there is no obvious difference between the two methods.

From Fig.2,it is evident that the ef ficiency of SH reconstruction is dependent on bothnmaxand the number of meshes to approximate irregular surfaces.Next,we will quantitatively discuss influences of the two factors via comparing traditional shape indices between SH-reconstructed shapes and actual ones.Considering that the target shapes are composed of 1500 sur ficial meshes,1280-face-based LBS particles are reconstructed withnmaxvarying from 1 to 15 by both Cartesian-and spherical-SH.In general,both Cartesian-and spherical-SH reconstructions have relative errors below 5%atnmax=5 forArandnmax≈10 for other indices,as shown in Fig.3.Therefore,a maximum reconstruction degree of 15 is suf ficient to well capture the main shape indices with relative error less than about 3%.

2.3.Mesh quality

The SH-reconstructed continuous surface is then discretised according toθandφ.The mesh starts with a regular icosahedron with its triangular elements divided into four identical triangles.The values ofθandφare determined by projecting the vertices on a sphere.By repeating this operation,we can discretise the surface into 20×4ntriangular meshes with the integern≥0.

Fig.3.Evolution of averaged relative errors of shape parameters for LBS particles with surfaces reconstructed with different maximum SH degrees for Cartesian-SH(a),and spherical-SH(b).

Fig.4.In fluence of surface discretization on particle shape parameters for the surfaces reconstructed with Cartesian-SH(a),and spherical-SH(b)(Weibull distribution plots with fitting lines).The vertical axis indicates ln(ln 1／P),and the horizontal axis indicates ln Sp,where P is the survival probability and Sp is the shape parameter.

Fig.5.Evaluation of surface mesh quality of LBS particles with the normalised element area,n A(a),and the length ratio,n H(b);and three kinds of surface meshes for a typical LBS sand(c).The length ratio is calculated as the shortest edge divided by its associated height.PDF denotes the probability density function.

To investigate in fluences of the number of sur ficial elements on approximating actual grains,we reconstruct 82 scanned LBS particles and 46 fragments with Cartesian-SH atnmax=15,since as indicated in Fig.2,it is shown that spherical-SH is not suitable for LBS fragments,many of which are not star-shaped.The reconstructed surfaces are discretised at different levels and then quanti fied with shape indices.The distributions of all shape parameters are very close for the surface meshes with 5120 and 1280 elements(Fig.4).Particle angular features are largely smoothed with 320 elements,which leads to larger roundness indices.Note that the quanti fication of roundness is all performed on the surfaces with 1280 elements by increasing or decreasing surface meshes.

The number of triangular elements and their shape quality have a great in fluence on simulation ef ficiency and precision(more details in Section 3).Two parameters are de fined to evaluate the mesh quality:the normalised element area,nA,and length ratio,nH.The surface meshes discretised from Cartesian-SH reconstructions have lower element quality with higher divergence for both the normalised element area and length ratio,simultaneously SHgenerated surface meshes have higher qualities than generalised marching cubes algorithm(Fig.5a and b).Remarkably,surface mesh sizes(Fig.5c)from marching cubes algorithm are visually less uniform than those from SH reconstructions,whereas spherical-SH generates better mesh qualities than Cartesian-SH,of which a vertex is shared by 6 triangles while by 4-8 triangles from marching cubes method.The lower mesh quality of Cartesian-SH is mainly caused by the surface parameterisation,mapping irregular surface onto a sphere,which leads to mesh distortions.

3.Single particle crushing simulations

3.1.FDEM

FDEM simulates an assembly of solid deformable bodies which can fragment(Munjiza,2004).FDEM models continuous phenomena (particle deformability),discontinuous phenomena(interaction and motion of particles),and the transition from continuum to discontinuum(particle fragmentation).We developed a subroutine to simulate the interface between elements in a commercial finite element method(FEM)solver Abaqus/Explicit 6.14(Dassault Systemes,2014).The cohesive interfacial element(CIE)degrades after a critical stress level is reached and eventually breaks with further energy input(see Appendix A).The contacts between discrete element objects are detected by general contact algorithm available in Abaqus.Individual sand particles are simulated with 4-node tetrahedrons(C3D4)bonded with zero-thickness and 6-node CIEs(COH3D6).The ef ficiency of FDEM in simulating single grain crushing has been benchmarked with an in situ CTmonitored experiment in a recent paper(Wei et al.,2019).

3.2.Single particle crushing tests

We model crushing behaviour of the LBS particle shown in Fig.1,the volume of which is altered to that of 1 mm-diameter sphere.The particle surface is reconstructed with both Cartesian-and spherical-SH atnmax=15.The reconstructed surface is discretised with 80,320,1280 and 5120 triangles,respectively,and converted into tetrahedral meshes.Then,we perform single particle crushing simulation on these six particles between two flat platens.The material properties,loading conditions,and solver precisions are identical.Notably,for comparison with existing XFEMstudy,all the material parameters are the same as those in Druckrey and Alshibli(2016)without any calibration.The simulated quartz sand has Young’s modulus of 94.4 GPa,Poisson’s ratio of 0.118;the maximum interface strengths are 25.3 MPa,12.6 MPa and 8.7 MPa for tension,first shear and second shear directions,corresponding to mode I,II and III fracture types of classical Grif fith fracture mechanics,respectively,for crack initiation.While the crack evolution is energy-based independent-mode and only mode I energy release rate(GI=100 N/m)is needed(Dassault Systemes,2014).The particle is placed with the smallest principal dimension perpendicular to the loading platens moving at a constant rate of 2 nm per step.More information about time step of FDEM can be referred to Munjiza(2004).

Fig.6.Load-displacement curves for FDEM-simulated single spherical particle crushing with different mesh sizes.The unit of mesh size is mm.

Firstly,for mesh size effect on breakage behaviour,the forcedisplacement curves of 3D spherical particles with various mesh sizes are provided in Fig.6.The particles generated with less surface meshes tend to have a higher initial stiffness(Figs.6 and 7).More sur ficial meshes result in finer characterisation of contact morphology,leading to higher local contact curvature and stress concentration.The Cartesian-SH-based meshes have large fluctuations before the peak force,which are mainly due to the low mesh quality,especially for 80 and 320 elements.According to McDowell(2002),the average failure force of one-dimensional(1D)compression of single LBS grain is 59.01 N,which is comparable to our results in Fig.7.

The von-Mises stress distribution shows a higher stress concentration at the contact area before peak stress for the particles with a denser surface mesh(Fig.8a).The fracture patterns also vary largely for six simulations(Fig.8b).For the simulation with the densest mesh,fractures initiate around the contacts,and then rapidly propagate towards two loading platens.Many fragments exist at the contacts due to shear stress.The coarser meshes have much simpler fracture patterns since fracture can only propagate through interface elements.

In Druckrey and Alshibli(2016),crushing processes of two CTscanned grains are simulated.The ratio of average mesh size to grain diameter is 0.0375(=0.03 mm/0.8 mm),which is nearly identical to that(0.037=0.037 mm/1 mm)of our 1280-meshcomposed surface.By comparing mesh sensitivity of spherical particles with the same material properties in both studies,it is evident that our simulation does better in convergencies of initial elastic response.In McDowell(2002),average forces of failure of experiments on natural quartz sands are 33.12 N,59.01 N,and 149.05 N for sizes of 0.5 mm,1 mm,and 2 mm,respectively.Surprisingly,the failure force of 0.8 mm sphere in Druckrey and Alshibli(2016)is about 110 N,which is signi ficantly higher.Meanwhile,the failure force of 1 mm sphere in our FDEM simulations is about 55 N,as shown in Fig.6.On the other hand,Wang and Coop(2016)experimentally studied 1D crushing behaviour of single quartz grains and concluded that different breakage modes existed,such as chipping,explosive,and the combined mode of the two.That is to say,real fractures in grain crushing are highly possible to cross with each other.Meanwhile,due to the characteristics of enriched elements of XFEM,one element can only have one propagated fracture,indicating that real fractures cannot be successfully simulated.Furthermore,in the simulations of Druckrey and Alshibli(2016),only one single flat fracture appears because of the dif ficulties in introducing density vibration in solid elements.Looking back the simulations in this study(Fig.8),FDEM simulations do produce more realistic fracture patterns.

Fig.7.Load-displacement curves for the single particle crushing tests of an LBS particle discretised in different surface element numbers for Cartesian-SH(a),and spherical-SH(b).

Fig.8.von-Mises stress distribution at the peak loading force(a),and the fracture patterns at 0.06 mm displacement(b)for an LBS particle discretised in different element numbers.

4.Conclusions

Simulating crushable irregular 3D particles remains a challenge because of the dif ficulties in reconstructing real particle surface and generating high-quality meshes.This study reconstructed smooth particle surfaces with SH in both Cartesian and spherical coordinate systems and discretised them into surface and volume meshes.We investigated the quality of generated meshes and their in fluence on FDEM simulation results for single particle crushing tests.The results lead to the following conclusions:

(1)Compared to other computational tools(e.g.DEM clusters and XFEM)to simulate single particle breakage behaviours,FDEM can not only effectively take particle shapes into consideration,but also generate more reasonable mechanical responsible and more realistic fracture patterns.

(2)For both Cartesian-and spherical-SH reconstructions for FDEM surface meshes,their mesh qualities are higher than that from commercially generalised marching cubes algorithm.The aspect ratios and roundness features are well captured at a maximum SH degree of 5 and 10,respectively,with relative error in shape indices less than about 5%.The spherical-SH reconstructions may generate unrealistic surface features,such as bumps and protrusions,especially for more irregular particles.However,Cartesian-SH may generate mesh distortion due to surface parametrisation.

(3)It is found that a higher mesh number results in a higher local contact curvature,leading to lower loading stiffness and higher peak strength for irregular particle crushing tests.Reducing the mesh elements may induce simple fragmentation patterns due to the less breakable interface elements.

Declarationofcompetinginterest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to in fluence the work reported in this paper.

Acknowledgments

This work was financially supported by Australian Research Council(Projects DP170102886)and The University of Sydney via SOAR(Sydney Research Accelerator)Fellowship.This research was undertaken with the assistance of the High-Performance Computer(HPC)service at The University of Sydney.

AppendixA.Supplementarydata

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2021.07.016.