1 Introduction

In industry, the purpose of applying DEM is commonly to gain a better understanding of how to optimise a manufacturing process, one which invariably includes dynamic granular flows. To be representative, the simulated particle numbers involved are typically many (tens or hundreds of) millions and the granular flow process of interest for simulation may have a duration of many seconds or minutes, e.g. sieving/screening, hopper/chute/belt conveying, crushing, milling, mixing, coating, tableting, etc. In such cases, even with massively parallelised computer hardware, accuracy of the particle geometry and hence particle-to-particle interaction often must be forfeited. Indeed, sometimes the governing physics is highly complex, especially for fine particles dominated by adhesive forces, so that non-spherical particle geometry is considered only of second-order importance and spherical shape approximation is adopted. Approximate matching of emergent bulk granular properties to experiments after lengthy application-driven calibration might be regarded as justification to proceed with further simulation to assist an industrial design process. An example might be multiphase fluid-coupled particulate flows in fluidised beds with free surface as discussed by Sakai [1]. To achieve run completion in practical timescales demanded by the engineer, simplifying assumptions are made to represent the particle shape and shape distribution and the size distribution of the granular medium to be modelled. This compromise is especially necessary when dealing with natural systems of rock fragment grains.

The packing process and packed topology of granular systems have been a target application of DEM methods ever since they were invented by Cundall and Strack [2]. The packing of particles and their packed characteristics (e.g. see German’s 1989 treatise on the subject [3]), as compared with granular flows in general, are a branch of granular problems that is more amenable to complex particle shape simulation. Improvements in computation power, parallel technologies and advanced, for example, contact algorithms in DEM and FDEM are, at last, sufficient to have impact on various packing applications used in industry. The need for and a specification for how to address the important role of non-spherical shape in gravity-dominated, for example, coarser particulates behaviour were highlighted a decade or so ago [4, 5]. For applications where the effects of complex shapes are of critical importance to fully represent a process with sufficient accuracy, extra computational effort must be expended on capturing the highly resolved shape; this requires comparatively smaller surface shape-defining discretised elements. Many industrial packing structures are physically generated by a relatively fast container-filling deposition process that brings coalescing moving particles rapidly in real time to a rest state. If in addition, they do not involve too many particles, they make ideal target problems because the associated smaller time steps required for the smaller elements and accurate contact dynamics can still be achieved within reasonable run times. As discussed later in the paper, this is even more so where the pack of complex-shaped bodies is itself a mono-shaped and mono-size system, for example, in catalytic reactors. For catalytic packed beds, the motivation relates to the need for an accurate capture of the topology of the void space for optimisation of gas flows with minimal and homogeneous pressure drop and to the contact heat transfer and potential fragility of the catalyst support grains/pellets. In the specialised case of coastal structures, the pack is built with concrete armour units and the stability of a one- or two-layer sloping pack under an oscillatory forcing wave–structure interaction is the problem. Both hydraulic and structural stability and flow through the pack become the targets of fluid–solid coupled simulation, aimed at improving breakwater armour layer system design.

In natural systems of packed rock fragment grains, as required for a fundamental study of fluid flow and soil or rockfill behaviour, progress has been made with highly sorted and rounded sand deposits. A numerical process DEM model using clustered sphere particle shape representation [6] of a sandstone’s structure was used to target fluid flow behaviour in a simple idealisation of a sandstone. The problem of assigning representative sizes and forms/angularities to the population of grains selected for the simulation domain, i.e. size distributions and shape distributions, surface texture, etc., for a specific geotechnical, for example, rockfill design application is still largely an outstanding one, especially for geomaterials containing angular grains. Even with an extensive shape library of numerically realistic rock fragment grain shapes [7], it is not trivial to determine, based on objective criteria, the proportions and sizes of each numerical grain shape to be deposited for a simulation designed to deliver an expected accuracy. This is a topic that requires further research. The effort to make the materials more representative is working in opposition to the extra computational expense as the application moves away from mono-size and mono-shaped niche problems. As an alternative to using DEM for sedimentation process models to create soil and rock fabrics and thereby investigate the effect of grain shape on geomaterial strength and flow properties, direct sampling of the actual grain fabric from rock cores, for example, by X-ray micro-CT is an option that is growing in popularity with digital imaging. This approach is discussed by Jia and Garboczi [8] in relation to shape acquisition and representation of realistic angular rock grains.

In this paper, we present a short summary of available numerical methods for modelling packing processes before focussing on two niche applications requiring high accuracy in the representation of complex shape, as made possible by FDEM.

2 Modelling methodologies for the effect of shape on packing

Packing algorithms for non-spherical particles were recently included in a general review that included granular material flows and some fluid–particle coupling schemes [9]. The modern algorithms that have superseded the once popular type of frozen-once-placed ‘ballistic deposition’ algorithms previously considered the most appropriate for spheres or for sphere composites [10] can be divided as follows:

  • Purely geometric These include random space filling and collective rearrangement, for example the ACS algorithm of Torquato and Jiao [11] used to obtain densest random packings of Platonic and Archimedean solids. Voxel-based methods [12] began as almost entirely geometric methods, the attraction being their ease for capturing any shape. Realism of such approaches has improved with the addition of some DEM-like mechanical interaction features to create ‘DigiDEM’ [13].

  • Geometrically determinate contact normal The DEM approach is based on extensions from the original spherical-based approach with geometrically determinate contact normal. Geometric forms such as ellipsoids and superquadrics [14] and as widely exploited in DEM technologies of Cleary’s CSIRO group, and termed superellipsoids [15], are invoked for the particle shapes, see also the review article [16].

  • Clustered spheres DEM based on clustered overlapping variable size spheres (clumped spheres) was made possible with geometry matching algorithms [17, 18]. These methods are the backbone of commercial simulators such as PFC3D (by Itasca) and EDEM (by DEM Solutions) and retain good computational speed performance at the expense of accuracy loss in capturing angular shapes. About ten spheres or less per complex particle geometry modelled has been suggested to be optimal for rock fragment processing problems. Simulation results of conveying and processing can seem very impressive in movie visualisations. However, in some cases, whereas the mechanical computation solution is based on clustered sphere interactions, the movie simulation summary of results often shows particle or object surfaces as a rendering of more angular highly resolved CAD-type particle/object surfaces. Examples may include conveying of bulk components like bottles or fastening brackets. Alternatively, shape library particles with faces and facets may be presented to look realistic. The rendered appearances are ‘draped over’ the clustered sphere geometry for which the DEM is solving. This normally hides the approximated bumpy surface that is in use and which reduces the accuracy of the sliding mechanics being simulated and their ability to represent sharp edge and corner interaction behaviour. Smooth joint approaches can be applied to clustered spheres to compensate for artificial bumps of the clustered sphere particle surfaces. DEM for ‘spherosimplices’ such as spherocylinders or spheroplates, for example, formed by moving a sphere along a line or perimeter can approximate certain shapes and still provide simple and fast contact interactions solutions [16].

  • Polyhedral DEM The particle shape is represented with surface discretisation where surface mesh element resolution can be set to capture the geometry considered necessary. Contact forces can be calculated between arbitrary-shaped convex and concave polyhedral particles with repulsive force formulations using ‘overlapping volume’ or ‘common plane area’ or ‘contact normal overlap depth’. Simple and fast solutions and linear scaling have been reported [19]. Voronoi tessellation methods have been used [20] to generate polyhedral shapes for ballast and rockfill applications, and a bounding box stage has been used to accelerate contact searches in some recent polyhedral DEM developments [21]. ROCKY is a commercial code developed from the polyhedral DEM approach that includes efficient breakage models.

  • FDEM (Combined FEM and DEM) The combination has been realised according to two main approaches to meshing. In the first, arbitrarily complex, for example, convex and concave particle/object shapes are constructed from combinations of primitive convex polyhedra. The bodies can be given rigid or deformable properties. An example is the 3DEC™ code by Itasca. Second, FDEM exploits surface triangular and volumetric tetrahedral meshes directly at the appropriate resolution to approximate the effects of any shape and these build from the pioneering suite of concepts and algorithms introduced by Antonio Munjiza. An important contribution to computing contact interactions of a generalised type is the application of a distributed ‘potential contact force’ algorithm [22]. It is possible to retain a fully deformable description of the multi-body system or alternatively to apply a rigid body approximation for each body [23], depending on whichever is best suited to the application and accuracy needed. The general-purpose codes ELFEN by Rockfield and IRAZU (specialising in geomechanics) by Geomechnica are perhaps the best known commercial FDEM codes and there exist several research groups and research laboratory FDEM codes, such as HOSS (Los Alamos National Laboratories) and open source Solidity (AMCG, Imperial College). Numerical models with similar capabilities to these FDEM codes have been created and are sometimes referred to as multi-particle FEM (MPFEM), especially in the field of powder metallurgy, where the focus is on plastic particulate properties [24]. Compared with FEM, the FDEM approach is more suited to handling fracture and fragmentation together with a higher number of multi-body contacts although traditional multi-purpose commercial FEM packages such as ABAQUS–Explicit and ANSYS (with RBD) are increasingly being augmented with contact detection and collision mechanics formulations, borrowing concepts from DEM [25] to create FDEM type capabilities. It would appear that a convergence is occurring as the FEM solids community augments their capabilities in recognition of the importance of FDEM.

For many applications with natural rock fragment shapes, the problem has been to introduce facetted angular particles and to solve the contact mechanics problem. Polyhedral DEM solvers handle better the angular contact interaction of crushed rock than clustered sphere representations. To represent an angular source material in polyhedral DEM, the BLOKS3D DEM program [26] creates typically ten or so different digital simple facetted particles by iteratively slicing separate grains until optimally matching predefined descriptors of angularity, flatness and elongation of each of the ten real aggregate particles, from which a representative digital sample is populated. This is an important development to encourage objective shape distribution representation.

A challenge for the FDEM methods, especially transient dynamic deformable simulations that track the stress inside the grains, is that the methods are relatively more expensive in terms of computational time, which limits the number of particles that can be considered or granular flow process time that can be modelled. Where the particles are known to have both smoothly curved and flat faces while also having sharp edges, such as cylinders, the number of tetrahedral elements to represent this particle geometry accurately often appears prohibitive. In this paper, illustrative results of the parallelised FDEM code Solidity [23] developed at Imperial College London by the Applied Modelling and Computational Group are presented in the context of two applications where close replication of particle geometry and resultant packed properties is essential for effective exploitation by the respective groups of industry clients.

3 Applications of FDEM using Solidity

Two very different applications featured in this paper are special cases where the system of interest can be meaningfully modelled with several hundred or a few thousand particles and the individual solid particle geometries themselves can be discretised using about 102–104 tetrahedral elements per particle.

3.1 Pellets for fixed bed catalytic reactors

Fixed bed reactor pellets, sometimes referred to as catalyst supports, will typically be made from porous ceramic materials to enhance reactions of the catalyst (Fig. 1) [27]. They will be shaped to offer a high surface area, to enhance reaction rates as pressurised fluids flow up and through the high temperature packed bed, often within steel cylindrical container tubes. However, for an ideal packed bed structure created from filling the cylinder, the random-like structure must not encourage localised channelised flows and not cause heterogeneous and local volumes with retarded gas flows, and therefore, a great deal of research into packing is in fact a precursor to the CFD modelling of heat and mass transfer through the pack [28]. In a development of rigid body dynamics with a hard body contact model for pellet packing, Moghaddam et al. [29, 30] and Garcia et al. [31] introduced background to packing and CFD pressure drop modelling in the industrial context of packed bed reactors during gas reforming. A widely reported difficulty encountered in CFD packed bed modelling is that the contact points between pellets lead to intangible meshing of the void space and lack of computational stability associated with non-robust meshing around contact points between pellets. Many strategies include shrinking particles slightly to create a gap, expanding the particles to give an overlap, expanding the contact to make a bridge or simply flattening contacts [32]. Moghaddem et al. [30] applied a ‘contact point treatment’ with volume inflation to overcome this problem and derive a computationally robust mesh for fluid flow modelling. This enabled them to use ANSYS-RBM workbench (with Fluent 16.2) software to compute fluid flow within the void topology bounded by the solid skeleton surface which was modelled with a rigid body dynamics (RBM) solver.

Fig. 1
figure 1

Steam methane reforming principles and catalyst pellet packs in reactor tubes [27]

Full size image

Catalyst pack geometry has also been generated by AMCG’s FDEM code Solidity, and a full coupling to a Navier–Stokes CFD solver, Fluidity, is achieved using an immersed body method capable of giving high velocity gas flow and pressure drop simulations for fixed bed problems. There are two distinct advantages of this approach. First, a high accuracy packing capability with Solidity is available for any complex concave or convex shapes (see Fig. 2) and a validation study by Farsi for spheres and trilobes [27] based on experiments showing good agreement is summarised in Figs. 3 and 4 and the ability to pack any shape illustrated in Fig. 5. Second, the ‘contact point problem’ does not arise as the adaptive fluids mesh resolves down to define the intricate topology [31]. Small but insignificant overlaps in the solids mesh always exist to regulate the contact forces through penalty function methods, and this has the benefit of preventing ill-formed volumetric meshes for fluid flow in a convenient way. A FEM mesh for fluidity (with adaptive mesh optimisation) completely separate to the solids mesh is generated for computing fluid pressure and transport. This mesh envelops the entire computational domain and tracks the solid volume fraction and the fluid flow. Solidity is coupled to Fluidity for both one-way and two-way coupled fluid structure interaction problems as explained in [33,34,35,36].

Fig. 2
figure 2

Packing structures created by the deposition of spheres and trilobes using Solidity [27]

Full size image
Fig. 3
figure 3

Packing density fluctuations in spheres: radial packing density distributions extrapolated from the X-Ray CT scans (in black) [48] and from the FDEM Solidity numerical simulations in a cylindrical container (red) [27]. (Color figure online)

Full size image
Fig. 4
figure 4

Packing density fluctuations in trilobes: radial packing density distributions extrapolated from the X-Ray CT scans (in black) from [48] and from the numerical simulations in a cylindrical container (red) [27]. (Color figure online)

Full size image
Fig. 5
figure 5

Packing structure at bottom of a cylindrical container simulated with FDEM in Solidity for ‘lettering pellets’ made from CAD file extrusion; to illustrate packing of a highly concave-shaped particle, and close-ups showing the tetrahedral mesh elements

Full size image

In steam methane reforming (Fig. 1), heat transfer in the flowing gas and amongst pellets from the side flame heaters to the tube’s interior pellets must also be efficient to drive the endothermic reactions, heat flow being another target of CFD/solids simulation. Furthermore, breakage of pellets into fragments causes reduced flow efficiency; uneven flow and excessively hot spots can result which in turn can adversely affect the steel container walls. Stresses additional to those from self-weight can build up from thermal stresses, and especially by shrinkage of the steel cylinder itself, if the catalyst process plant is shut down and reverts from ~ 900 °C to ambient temperatures [27]. To better understand the role of pellet shapes and the strength and physical properties of the support medium on the intricate coupled behaviour, the ability to model thermal conduction, shrinkage and possible (thermally assisted) fracture behaviour of the pellet pack would be a distinct advantage in a versatile tool for modelling packed bed reactors. FDEM can capture the size distribution of fragmented catalyst pellets. Ideally, the effect of any breakages, fines production and clogging could then also be evaluated with the coupled high Reynolds number fluid flows of the gas captured, for example, by an adaptive CFD code. Models for heat flow within pellets with arbitrary surface boundary conditions and across multi-body contacts have been developed for FDEM Solidity (see Fig. 6) [37]. FDEM models in Solidity are also able to simulate thermal stresses related to volumetric shrinkage and expansion strains within the pellets leading to thermal fracture [38, 39]. Turning to stress transmission in general, Fig. 7 shows stress chains resulting from a gentle compression of packed cylinders. For a simpler 2D illustration of the role of stress chains in FDEM simulations, see Fig. 8 [40].

Fig. 6
figure 6

Heat conduction from a hot sphere (on right) contacting a cold sphere on left, model implemented in FDEM Solidity, by Joulin [39]

Full size image
Fig. 7
figure 7

Vertical slice and 3D view showing differential stress computed inside solid cylindrical pellets packed in a cylinder with FDEM Solidity. Stress chains and vulnerability to crushing or tensile fracture can be examined depending on stress components inspected, [27]

Full size image
Fig. 8
figure 8

Packing and stress chain structure of high aspect ratio ellipses with friction coefficient of 0.5 in 2D FDEM, see [40] for detailed discussion

Full size image

For simulations that capture internal particle deformation and stresses, computational demand is about an order of magnitude greater than with rigid particle approximations. However, it is possible to import packed geometry meshes created using rigid modelling and to bring them to equilibrium for later loading in a fully deformable FDEM model. This makes for an efficient work flow to examine stresses within larger realistic pack systems (Figs. 7, 8). A further advantage of FDEM is the potential to apply different constitutive models for the deformation response to stresses. Solidity supports validated 2D and 3D fracture models [27, 41, 42] applicable to granular breakage and fragmentation in multi-body loading (Figs. 9, 10).

Fig. 9
figure 9

3D fracture modelling [42] and 2D fracture models of disc with holes that improve surface area but make pellets more prone to fracture [27] in FDEM Solidity

Full size image
Fig. 10
figure 10

Pendulum collision test with concrete dolos unit used in validation study of 3D fracture model left: velocity, centre: dynamic tensile stress, right: fracture visualisation [41]

Full size image

Alternatively, if the particle materials are elastic–plastic, a range of plasticity models have been implemented in FDEM Solidity [43] to predict ductile multi-body compaction or extrusion (Fig. 11) and local crushing in concrete. The large strain plasticity capability in Solidity is also applicable to catalyst pelletisation (and other powder compaction processes) as pellets are formed from compacted soft aggregates.

Fig. 11
figure 11

High-speed collision of two plastic spherical bodies [43, 49] showing large strain plasticity in FDEM Solidity and the capability to model compaction in tablet and pellet forming

Full size image

3.2 Concrete armour units

Harbour breakwaters are often constructed with concrete armour unit layers. When designed appropriately, these gigantic granular systems can withstand the wave forces for a predicted sea state or design storm. Breakwater armour system design is largely empirical and would be transformed by a better fundamental understanding of loads and restraining forces. Numerical models based on FDEM are leading the solids modelling effort as they can handle the complex shapes and contact forces most realistically. The units are very massive (10–100 t) and are cast with unreinforced concrete. The units available have widely varying shapes, and if left loose and poorly interlocked may become dislodged or simply rock under wave action. As concrete is weak in tension and may crack, the interlocking and unit robustness is of great significance. Most unit types are placed with irregular orientations by crane in rows and columns to a GPS system grid to yield a specified packing density.

Numerical models allow statistical parameters such as coordination number, contact forces and tensile stresses of the static pack to be analysed [44] (Fig. 12a–c). Random wave sea states have been modelled to interact with the FDEM breakwater topologies by applying the wave simulation code IHFOAM [45] to obtain the time history of velocities and pressures of the turbulent flows. After calibration, the forcing hydraulic lift and drag forces are obtained. The hydraulic forces are summed together with the solids forces governing motion in the FDEM code Solidity in an application termed the ‘wave proxy’ which has been used to predict stability and armour unit movements, for different armour pack characteristics, for 1 h of a 100-year return period storm [45, 46] (Fig. 12d, e). Innovative shapes of units are frequently proposed. To encourage their take up by designers and clients when there is no record of site performance, demonstration by scientific principles embodied in mechanical simulations and analysis tools may prove more effective than empirical model testing in a laboratory where scale effects may not be accounted for and detailed use of sensors and analysis can be expensive and time-consuming. For example, the design of units with specially added geometric features such as the sacrificial friction studs of Accropode II™ can be interrogated by advanced modelling (Fig. 12f, g).

Fig. 12
figure 12

Application of FDEM to breakwater stability prediction [44, 45]: a breakwater with Core-Loc units in row–column pattern, b model of breakwater section ~ 106 elements with statistics of maximum contact force, c tensile stress details, d pressure and velocity time history for wave–structure interaction sequence modelled with IHFOAM, e movement of armour units modelled with FDEM Solidity wave proxy, f Accropode II™ unit, g position of armour units before (white) and after (colours) in FDEM model of storm wave sequence. (Color figure online)

Full size image

4 Discussion and concluding remarks

To illustrate how simple it is to explore ideas of science and engineering with FDEM codes such as Solidity, two armour unit types were modelled as catalyst pellet size particles of equal volume and deposited from an identical evenly spaced array into a cylinder to compare random packing behaviour for two armour unit shapes used widely in coastal structures, Core-Loc and X-bloc [47]. When deposited in bulk with random deposition, the packing density and topology appear quite different for each unit, for identical friction coefficients; clearly there are fundamentally different packing phenomena affecting the packing characteristics for each unit (Fig. 13), but the collision dynamics in such small particle models involve only very small masses. Many scaling issues are of potential interest: at full scale, the dynamic stresses generated would be highly complex and damaging to the units. If captured with the deformable solver, the consequences for breakages and concrete fragments could be evaluated, provided sufficient computer resources could track the 3D fracture modelled response. When full-scale units are placed gently one by one in breakwaters using regular grid spacing and semi-random orientations, single layers made of both Core-loc and X-bloc types of units achieve similar packing densities and are found to perform with very similar stability under wave action. This suggests that the slinging conventions and placing by the crane operator have overriding importance in what are in fact highly non-random final positions, positions that are critical for these units to work and interlock effectively. The packing and interlocking properties from simulations from ‘random settlement bulk packing’ construction conditions such as those in Fig. 13 make for interesting comparison with simulations from contractor one-by-one and row-by-row grid-based construction. With FDEM solvers, a whole suite of analysis tools to interrogate old and new shaped armour unit systems for hydraulic and structural stability is in prospect.

Fig. 13
figure 13

FDEM Solidity simulation showing packing structure of concave particle geometries, Core-Loc™ left and X-bloc® right, widely used for single-layer armouring of breakwaters. Random orientation and spacing structure for particle entry of ~ 3 mm width concave-shaped particles, i.e. ~ 1000 times smaller than typical breakwater prototype sizes. Note: higher bulk porosity of X-blocs

Full size image

In summing up, it is proposed that we can now identify a distinct class of ‘FDEM-suited’ problems of complex-shaped particulates. The class of problems has two properties: firstly, a high accuracy in shape representation and contact mechanics is essential to the simulation user, and secondly, the number of particles is relatively few such that industrially useful simulation results are achievable in practical runtimes on affordable hardware—i.e. a company can improve the profitability of its industry/organisation by adopting such simulation tools. DEM has made its mark in particle technology industries already, but conventional DEM is likely to struggle in this FDEM-suited utility class of particulate problems. The applications illustrated here are catalyst packing behaviour and armour unit pack stability, and the most suitable software, we argue, will be of the FDEM type. For both these illustrated problems, packing characteristics can be indicators of structural stability and/or resistance to flows through the packs which can be readily investigated with sensitivity analyses. Computational power of affordable multi-core desktops is today sufficient for FDEM codes (with parallelised architecture and advanced interaction and search algorithms together with algorithmic efficiencies) to tackle industrial problems with computer runs of about 1 day or less for this shape-critical class of problem. With today’s technology, these FDEM-suited problems identified to date are: the relatively rapid packing of relatively few particles of complex shape that can be captured within the computational domain with, say, ~ 105–107 tetrahedral elements. Both are examples of what was termed ‘proactive use’ of shape information by Jia and Garboczi [8]. It is suggested that today, several industrial packing problems are already in this FDEM-suited utility class or will be in the very near future, including the two applications highlighted in this paper.