Constriction size distributions of granular filters : a numerical study

The retention capability of granular filters is controlled by the narrow constrictions connecting the voids within the filter. The theoretical justification for empirical filter rules used in practice includes consideration of an idealised soil fabric in which constrictions form between co-planar combinations of spherical filter particles. This idealised fabric has not been confirmed by experimental or numerical observations of real constrictions. This paper reports the results of direct, particle-scale measurement of the constriction size distribution (CSD) within virtual samples of granular filters created using the discrete-element method (DEM). A previously proposed analytical method that predicts the full CSD using inscribed circles to estimate constriction sizes is found to poorly predict the CSD for widely graded filters due to an over-idealisation of the soil fabric. The DEM data generated are used to explore quantitatively the influence of the coefficient of uniformity, particle size distribution and relative density of the filter on the CSD. For a given relative density CSDs form a narrow band of similarly shaped curves when normalised by characteristic filter diameters. This lends support to the practical use of characteristic diameters to assess filter retention capability.


INTRODUCTION
Granular filters are placed in zoned embankment dams and flood embankments in order to protect adjacent base soils which are vulnerable to internal erosion, in particular the watertight core material.Granular filters function by retaining base material within the narrowest constrictions in the void network between the filter particles.Particle retention within granular materials is also important in other areas such as water treatment (Yuan et al., 2012), bacterial filtration (Martin et al., 1996), sand production in wells (Penberthy & Shaughnessy, 1992) and leachate collection systems (Rowe, 2005).
Design criteria for granular filters typically take the form of a comparison between a base particle diameter and a filter particle diameter, where the filter particle diameter implicitly characterises the constrictions within the filter.For the protection of cohesionless base soils, Sherard & Dunnigan (1989) and ICOLD (2015) recommend Terzaghi's classic filter rule that D 15F , 4 D 85B , where D XF is the filter particle diameter for which X% of the material by mass is finer and D XB is the base particle diameter for which X% of the material by mass is finer.In Terzaghi's original patent application the rule is justified on the grounds that 'the pore size of a broadly graded filter comprises at maximum 1/5th of the diameter of the biggest grain of the finest fraction of the filter materials' (Fannin, 2008).This implies that if D 15F , 5 D 85B , then the D 85B will not pass through the pores of the filter (the design rule includes a safety factor).The same basis has also been used to justify the use of characteristic diameters in other areas of engineering, for example the internal stability of soils (e.g.Kézdi, 1979), permeability (Chapuis, 2012) and leachate collection systems (Yu & Rowe, 2012).Kenney et al. (1985) defined the 'controlling constriction size', D* c , of a granular filter as the diameter of the smallest constriction which a base particle is likely to encounter on a given flow path.They considered this to be the governing property of a granular filter.Kenney et al. (1985) carried out an extensive series of filter tests using uniform cohesionless base materials and linearly graded filters of varying thickness with coefficients of uniformity, Cu = D 60F /D 10F = 1•2 to 12. Extreme hydraulic conditions and light vibration were applied to ensure that the largest possible base particles were transported through the filter.The largest particle which was transported through each filter was taken to be equal to D* c .Kenney et al. (1985) related D* c to the characteristic diameters D 0F , D 5F and D 15F .For the most uniform filters with Cu = 1•2, they found D* c % 0•18D 0F , whereas for filters with Cu ! 3, they found D* c % 0•25 D 0F .When normalised by D 5F and D 15F it was found that D* c was similar for all filters regardless of Cu value, with upper limits to the values given by D* c 0•25 D 5F and D* c 0•20 D 15F .The results using D 15F were found to give slightly more scatter than D 5F , but were still sufficiently good to recommend the continued use of Terzaghi's D 15F /D 85B 4 ratio for design for filters with Cu 6. Limited experimental work has been carried out for filters with Cu .6, notable exceptions being Sherard et al. (1984) and Lafleur et al. (1989), both of whom found that D 15F remained a suitable characteristic filter particle diameter.However, based on permeability data, Lafleur & Tétreault (1986) found that relative density was more important for filters with Cu .6 than for more uniformly graded filters, suggesting that more research is needed in this area.
Analytical approaches to filter design in geomechanics have focused on the 'inscribed circle method', which is used to calculate the full constriction size distribution (CSD) and to relate this to the retention ability of the filter (e.g.Silveira et al., 1975;Indraratna et al., 2007).The CSD is the cumulative distribution of constriction diameters (by number) within a filter.The inscribed circle method pioneered by Silveira et al. (1975) calculates the CSD by considering the largest 'constriction' circle which can be inscribed in a highly idealised two-dimensional (planar) representation where both the filter particles and the constrictions are represented as circles, as shown in Figs 1(a) and 1(b).The 'dense' state constriction diameter, D cD , is given by the largest circle which can be inscribed between three mutually touching filter particles (Fig. 1(a)).The 'loose' state constriction diameter, D cL , is given by the equivalent diameter of the maximum area, S c,max , between four touching filter particles, as shown in Fig. 1(b).The probability of a filter particle forming a constriction (i.e.acting as a vertex in a combination of contacting spheres) is related to the probability of its occurrence in the particle size distribution (PSD) (either by number, surface area or mass).Silveira et al. (1975) related the probability of a particle of a given diameter forming a constriction to the PSD by number.The Silveira et al. (1975) method was further developed by Locke et al. (2001), who proposed that the CSD should be related to the PSD by surface area and should vary linearly with relative density according to the relationship where P is the fraction smaller of the CSD (i.e. the current point of the curve on the y-axis P = 0 at the smallest constriction and P = 1 at the largest constriction); D cP is the constriction size at P; R d is the relative density; and D cD , D cL are the dense (three-particle) and loose (four-particle) constriction diameters for P, respectively.Locke et al. (2001) also proposed using a regraded PSD to calculate the CSD so that loose fines, which do not contribute to the constriction-forming skeleton, and coarse particles, which are 'enmeshed in a matrix of fines', are not considered.These coarse particles are removed from the PSD as any constrictions between them will be filled with finer particles, so will not be true constrictions.Kenney et al. (1985) used the Silveira et al. (1975) inscribed circle method in their study into the controlling constriction size for granular filters.They used a one-dimensional infiltration model of a void network to determine the controlling constriction size (i.e. the maximum size particle that can be transported through a filter of a given thickness).For filters of at least ten layers with the same D 0F and with coefficients of uniformity ranging from Cu = 2 to 12, Kenney et al. (1985) found that the controlling constrictions were similar for all filters, in agreement with their experimental work.Noting that the method used takes the PSD by number to calculate the CSD, they attributed the similarity to the relatively large number of fine particles dominating the CSD and therefore the controlling constrictions.The analytical and experimental work of Kenney et al. (1985) formed the basis for the design rule developed by Kenney & Lau (1985).Locke et al. (2001) extended the work of Kenney et al. (1985) by replacing the one-dimensional model of Kenney et al. (1985) with a three-dimensional cubic model allowing greater degrees of freedom for the permeating fines, as proposed by Schuler (1996).Based on this, Indraratna et al. (2007) proposed that the controlling constriction size, D* c , should be equal to D c35 , the 35th percentile from the Locke CSD calculation method.Indraratna et al. (2007) argued that the advantage of calculating D* c from the CSD, as opposed to taking a characteristic value from the filter PSD, is that the complete PSD shape and relative density are taken into account.Overall, this method predicts that, as filters become more widely graded, the CSD will also become more widely graded.Note that the regrading rule proposed by Locke et al. (2001) is apparently not applied in the Indraratna et al. (2007) method.
Typical permeameter data can only reliably give values for D* c and cannot be used to define a CSD.The above analytical approach considers idealised combinations of spherical particles and hitherto the validity of the underlying assumptions has not been directly assessed.Advances in numerical modelling (discrete-element method (DEM)) and image analysis (micro-computed tomography) mean that constriction topology can be examined at the particle scale and a CSD can be determined (Ketcham & Carlson, 2001;O'Sullivan, 2011).The resultant data enable a scientific re-examination of the geometrical hypotheses that date from Terzaghi's work and underlie modern filter design (e.g.ICOLD, 2015).
Here the DEM is used to generate numerical samples for which CSDs can be determined to enable a re-examination of the conclusions developed from the inscribed circle assumption.The algorithm proposed by Reboul (2008) to generate CSDs is introduced.The inscribed circle  (Silveira et al., 1975); (b) idealised two-dimensional, four-particle arrangement used in the inscribed circle method (Silveira et al., 1975); (c) to (f) images of weighted Delaunay constrictions inscribed between three-dimensional DEM particles assumption is then critically analysed by comparing CSDs obtained from the method proposed by Locke et al. (2001) and Indraratna et al. (2007) with CSDs obtained using Reboul's algorithm.The variation in the shape of the CSD with PSD and density is then considered.Finally, the idea of using a characteristic diameter as a means to relate particle sizes and constriction sizes is quantitatively assessed.

ANALYSIS APPROACH
The DEM is a numerical approach that was developed by Cundall & Strack (1979) to enable particle-scale simulation of the mechanical behaviour of granular materials.In a DEM model each particle and each particle-particle contact is explicitly considered.Ideal particle geometries (spheres) are used to control the computational cost of the DEM simulations.The effect of particle shape on filtration is unclear and it is typically not considered in design criteria, although Wu et al. (2012) found that filters comprising glass beads gave a similar experimentally derived CSD to rounded sands, whereas angular sands gave a wider range of constriction sizes.The data directly generated from a DEM simulation include the particle positions, the contact locations and the contact forces.Shire et al. (2014) showed that these data can be used to assess the internal stability of granular filters.The current study extends this earlier contribution by using DEM data to quantify constriction sizes.
The DEM simulations were carried out on cubic samples using a modified version of the open-source DEM code Granular LAMMPS (Plimpton, 1995).Periodic boundary conditions were used to create a sample which is effectively of infinite size and is free from the boundary effects associated with rigid boundaries, allowing particle numbers to be kept to a reasonable level.A Hertz-Mindlin contact model was used and the simulation input parameters were Poisson ratio, ν = 0•3, shear modulus, G = 27•0 GPa, and particle density, ρ = 2670 kg/m 3 , which are approximately equal to experimentally derived values for spherical glass beads used by Barreto (2009).Particles are initially placed in random, non-touching positions within the periodic cell using an in-house placement code, creating a homogeneous, high-porosity sample, and then applying isotropic compression with a uniform strain field following the algorithm proposed by Cundall (1988).The target isotropic stress level was p′ = (σ′ 1 + σ′ 2 + σ′ 3 )/3 = 50 kPa where σ′ 1 , σ′ 2 , σ′ 3 are the three principal stresses.Samples created in this way are homogeneous and isotropic.Simulations were terminated when the mean normal stress reached the target level and the coordination number (the number of contacts per particle) remained constant for 20 000 simulation cycles.
Referring to Table 1 and Fig. 2, 18 of the DEM samples created had linearly graded PSDs with Cu = 1•2, 1•5, 2, 3, 4•5  or 6, a further ten of the samples had bilinear PSDs with Cu = 1•5 to 2•6.Each of the PSDs is presented normalised by the smallest filter particle diameter (D 0F ).The PSDs of each of the bilinear samples follow a linear PSD with Cu = 3 to a diameter D X , where X is specified by the number in the sample name (D 5 for BL5, D 10 for BL10, and so on), and a PSD of Cu = 1•5 thereafter.Up to three coefficients of friction (μ) were used for each grading considered (i.e.μ = 0•0, 0•1 and 0•3) to give a range of void ratios, and the resultant samples are termed 'dense', 'medium' and 'loose', respectively.μ = 0•3 was selected for the loose samples as this is approximately equal to the value reported for physical glass beads by Barreto (2009).The number of particles considered increased with Cu and varied between 8262 and 59 183.As detailed by Shire (2014), DEM simulations with PSDs with Cu = 1•5, 3 and 6 were repeated using samples with fewer particles to confirm that a representative element volume (REV) had been achieved.Referring to Table 1 it is clear that for a given μ value, the void ratio at p′ = 50 kPa falls as Cu increases.The maximum coefficient of uniformity considered was Cu = 6 to limit the computational time required for simulations, as discussed above.The approach used here to determine the constriction sizes and hence the CSDs was the weighted Delaunay tessellationbased algorithm originally proposed by Al-Raoush et al. (2003) and Reboul (2008).A two-dimensional schematic diagram of the algorithm is included as Fig. 3.A threedimensional weighted Delaunay tessellation was formed with the tetrahedra vertices being located at the particle centroids; during the creation of the tessellation the vertices were weighted by their corresponding particle radii (Edelsbrunner & Shah, 1996).The weighted triangulation was achieved using the CGAL Regular Delaunay algorithm (CGAL, 2013).Referring to Fig. 3(b), inter-void constrictions are identified using the faces of the tetrahedra.On each face the constriction diameter is taken to be the diameter of the smallest circle that can be inscribed between the three particles that form that face or between the particles forming the face and any other non-vertex particles which cross over the tetrahedron face (Reboul et al., 2010).
In order to avoid over-segmentation of the void space, adjacent tetrahedra were merged using a criterion proposed by Al-Raoush et al. (2003) and Reboul et al. (2008).This can be explained by reference to a simple two-dimensional analogy of uniform particles on a square grid, as illustrated in Fig. 3. Fig. 3(a) illustrates the constriction locations.
The weighted Delaunay triangulation to identify the voids is presented in Fig. 3(b).As illustrated in Fig. 3(c) if constrictions are located at each of the triangulation edges (equivalent to tetrahedron faces in three dimensions), the voids are essentially over-segmented, that is, a space that would intuitively be considered a single void is subdivided.As illustrated in Fig. 3(d) circles (spheres in three dimensions) that are tangential to all the particles forming a Delaunay cell are identified.Then the overlap between this tangent sphere and each of the tangent spheres in adjacent Delaunay cells is calculated and if this overlap exceeds a user-specified value, the cells are merged to form a single void cell, as illustrated in Fig. 3(e).
The choice of user-defined merging overlap value is subjective; a smaller value results in the formation of fewer, larger voids.Although this subjectivity is indeed a limitation of the algorithm, it is important to realise that as soils have a continuous network of interconnected voids the exact definition of boundaries between voids in the system is itself inherently subjective.The variation in the calculated CSD with the merging overlap selected can be appreciated by reference to PSD.For both Cu values, as the merging overlap increases fewer cells merge and more constrictions are measured; the constriction sizes also increase with merging overlap.The effect of changing the merging overlap is less pronounced for the smaller than the larger constrictions.Typically a merging overlap of 50% has been used in the analyses presented here; however, where appropriate more than one merging overlap value has been considered to account for the subjectivity (2010) compared CSDs generated with the triangulation approach to those calculated using the inscribed circle method.They considered loose and dense DEM samples of soils with uniformly to moderately graded PSDs (Cu = 1•67 and 3•75) and they also found that for both filter gradings that the inscribed circle CSDs for R d = 100% gave poor agreement and were much finer than the weighted Delaunay CSDs for their dense DEM samples (μ = 0).They did find that the inscribed circle method for R d = 0% gave a reasonable match to the loose DEM samples, and the deviation from the results here is probably because they used a higher inter-particle friction coefficient for their loose samples (μ = 0•7).In addition to the study of filtration, the weighted Delaunay method has also been used for the prediction of permeability in porous media (Gao et al., 2012).

Evaluation of idealised constriction configurations
The data that can be generated using this approach enable a re-examination of the hypothetical constriction configurations considered in prior studies and illustrated in Figs 1(a) and 1(b).These assumptions are embedded in the inscribed circle method proposed by Silveira et al. (1975) and Locke et al. (2001).Both the inscribed circle method and the weighted Delaunay method give data on the CSD as well as the particles bounding each constriction.Direct comparison of data generated using these approaches is therefore well suited to a re-evaluation of the use of idealised configurations when generating a CSD.The inscribed circle method was implemented using Matlab as described by Shire (2014).The implementation was successfully validated against an example from Raut (2006).
Figure 5 compares the CSDs obtained using the inscribed circle method with those obtained using the weighted Delaunay method for Cu = 1•5 (Fig. 5(a)) and Cu = 6 (Fig. 5(b)).The weighted Delaunay method clearly generates a much narrower range of CSDs than the inscribed circle method in both cases.For the uniformly graded material with Cu = 1•5 shown in Fig. 5(a), the weighted Delaunay CSDs plot approximately between the inscribed circle CSDs with R d = 30% and R d = 60%.Within these limits, the shape of the CSDs is in agreement.The densest possible arrangement of the inscribed circle method (R d = 100%) which assumes all constriction configurations are three-particle (Fig. 1(a)) gives constriction diameters that are much smaller than those obtained using the weighted Delaunay method for the 'dense' sample.The 'dense' weighted Delaunay CSD was created considering frictionless, perfectly spherical particles and so it represents an upper bound that is denser than any sample that could be physically created.Thus, it can be concluded that a fabric with constrictions consisting entirely of three mutually touching spheres is unlikely to occur in reality.Also, it is unclear whether it is possible for spheres to form a stable fabric that would be loose enough to match the loosest possible arrangement (R d = 0%), as the 'loose' DEM sample was created using an interparticle friction of μ = 0•3, which prior research suggests is close to the upper limit of mineral friction, as discussed by Huang et al. (2014).Huang et al. (2014) also showed that using a value of μ !0•5 can give a material response that is not representative of soil behaviour.The R d = 0% data for the inscribed circle method take the constriction formed by the four-particle configuration using the largest particles (Fig. 1(b)) to be D c100 ; the data presented here indicate that that assumption significantly overestimates the constriction sizes.
The narrower range of CSDs obtained using the weighted Delaunay method may be due in part to the use of perfect spheres, which tend to give a smaller void ratio range than real soils.However, even taking this into account, the weighted Delaunay result does not support the use of equation ( 1) to model the variation of the CSD with relative density.Figs 1(c)-1(f) illustrate representative particle-constriction configurations obtained in the DEM models.These are presented adjacent to the assumed three-particle and fourparticle configurations (Figs 1(a) and 1(b)) to illustrate that, although these idealised configurations occur in some cases, for example Figs 1(c) and 1(d), other constriction configurations do not resemble these idealised combinations, for example Figs 1(e) and 1(f).
Figure 5 (b) shows that the agreement between the weighted Delaunay and inscribed circle methods reduces as the Cu increases to Cu = 6.Although both methods give the same D c0 % 0•155D 0F , the inscribed circle method results in a more widely graded CSD than the weighted Delaunay method.The divergence between the weighted Delaunay method and the inscribed circle method for intermediate relative densities (R d = 50% and 70%) occurs for the largest 50-60% of constrictions in the CSD, where the inscribed circle method results in much larger diameter constrictions.The densest and loosest CSDs from the inscribed circle method diverge greatly from the weighted Delaunay result.The densest CSD generated using the inscribed circle method includes constrictions that are larger than those from any of the weighted Delaunay CSDs.
This reduction in the agreement between the two methods as Cu increases can be explained by examining the soil fabric.Fig. 6 gives the cumulative distributions of the largest, intermediate and smallest particles bounding constrictions for the inscribed circle and weighted Delaunay methods for the material with Cu = 1•5.The inscribed circle distribution is based on the assumption that the probability is proportional to a particle's surface area (Locke et al., 2001), whereas the weighted Delaunay distribution is directly measured from the numerical model.For this uniform filter, the surface area assumption is reasonably accurate, as shown in Fig. 6(a).The agreement decreases with increasing Cu and the cumulative distributions for the filter with Cu = 6 are presented in Fig. 7.It is clear that the particles assumed to form constrictions according to the inscribed circle method are consistently coarser than those from the weighted Delaunay.These data show that the assumption that surface area is a measure of how likely a particle is to form a constriction does not hold for the more broadly graded filters, as the coarse particles are separated by finer particles, as shown schematically in Fig. 7(b).This shortcoming was recognised by Locke et al. (2001), who suggested regrading well-graded filters by removing particles with D nF , D 2nF /4, where D nF is the filter diameter for n% smaller and D 2nF is the filter diameter for 2n% smaller.For linear gradings this equates to regrading the PSD to Cu = 4.It can be seen in Fig. 5(b) that this gives only marginally better agreement.
These comparison data indicate that the idealised configurations in Figs 1(a) and 1(b) are not representative of constriction topologies.This lack of agreement motivated a systematic re-examination of the influence of Cu and density on constriction sizes as well as a re-assessment of the characteristic diameters used in common design rules.

EFFECT OF COEFFICIENT OF UNIFORMITY AND DENSITY ON CSD
Figure 8(a) shows the CSDs obtained using the weighted Delaunay method with merging overlap = 50% for the loose, linearly graded samples, normalised by D 0F .All the CSDs for the linear gradings considered have generally similar shapes; referring to Table 1 and considering all packing densities, the CSD Cu varies between 1•48 and 1•64.As the filter grading increases from Cu = 1•2 to Cu = 3, the CSDs coarsen.However, for samples with Cu ! 3 the CSDs are very similar, with those for Cu = 4•5 and Cu = 6 being virtually indistinguishable, as shown by comparing the D c50 and Cu CSD values in Table 1.The simulations carried out with the bilinear gradings are also summarised in Table 1.Fig. 8(b) shows the CSDs for the bilinear PSDs given in Fig. 2(b).As with the CSDs from linear PSDs, the CSDs for the bilinear PSDs maintain a relatively similar shape.For BL5, BL10 and BL15 the CSDs lie between the limits set by Cu = 1•5 and Cu = 3, the CSDs becoming more similar to Cu = 3 as the volume fraction with Cu = 3 increases.However, as the proportion of Cu = 3 material moves to 25% and above the CSDs become very similar to, but slightly coarser than, the CSD of the linear filter with Cu = 3.The trends observed for the linear and bilinear gradings were also observed in additional analyses using merging overlaps of 0% and 100% and for the other density values considered.Relative density is an important factor in filter effectiveness (Lafleur & Tétreault, 1986), although it is not considered by empirical rules such as Terzaghi (Fannin, 2008) or Sherard & Dunnigan (1989).The effect of density is illustrated by comparing the loose, medium and dense CSDs for samples with Cu = 1•5 and Cu = 6 in Fig. 9.As shown in Table 2, the smallest constrictions are the same for loose and dense samples (% 0•155D 0F ), and the remaining fractile diameters (e.g.D c25 , D c50 , and so on) are consistently approximately 10% smaller in the dense sample than the loose, indicating that while the constriction diameters show some sensitivity to density, the shape of the CSD curve remains similar.
Comparison with characteristic diameters proposed by Kenney et al. (1985) Kenney et al. (1985) proposed normalising the CSDs by D 5F and D 15F to relate the controlling constriction size to the PSD.Fig. 10(a) illustrates the CSDs normalised by D 5F of the corresponding filter PSDs and it is clear that the CSDs from the samples with Cu ! 2 form a narrow band of curves, with the CSDs from the samples with Cu = 1•2 and 1•5 lying to the left of this (i.e.having smaller normalised constrictions).The controlling constriction size D* c % 0•25D 5F proposed by Kenney et al. (1985) is included on Fig. 10(a).Although the CSD alone does not allow the controlling constriction size to be calculated, qualitatively it can be considered that similar normalised CSDs would translate into similar values of D* c for these filters.Fig. 10(b) considers the CSDs normalised by D 15F and in this case the samples with Cu 3 form a narrow band, with more widely graded filters having progressively finer normalised CSDs.The controlling constriction size D* c % 0•2D 15F proposed by Kenney et al. (1985) is also marked on Fig. 10(b).Fig. 11 shows the CSDs normalised by D 15F for the bilinear PSDs.As was the case with the linear PSDs with Cu 3, the normalised bilinear filter CSDs maintain a similar shape and form a narrow band of curves.
Figure 12 The fact that the median constriction diameter depends on the chosen merging overlap highlights that the definition of the actual CSD is ambiguous.However, although a unique CSD cannot be defined, there are similar trends in CSD variation with Cu irrespective of the chosen merging overlap used to define the CSD.This lends support to the approach of Kenney et al. (1985), who proposed using an experimentally determined proportion of a characteristic filter diameter to obtain D* c , rather than attempting to explicitly calculate the full CSD for design.However, the sensitivity of the CSD to relative density shows that the PSD alone is insufficient to describe the controlling constriction at all density levels.
A comparison between the normalised median constriction sizes for each normalising parameter (D 0F , D 5F and D 15F ) is given in Table 3.This confirms that D 5F gives the narrowest band of normalised values for more widely graded filters with Cu ! 2, whereas D 15F gives the narrowest band of normalised values for more uniform filters with Cu 3.
Both have a similar overall scatter for the samples analysed, measured by the coefficient of variation (standard deviation divided by mean).For linearly graded materials with Cu 3 it was noted above that D 15F was the best normalising parameter.As shown in Table 3, the same pattern is found for the bilinearly graded materials, with the CSDs normalised by D 15F giving a narrower band of CSDs than the normalised CSDs of either D 5F or D 0F .

CONCLUSIONS
Inter-void constrictions are a key characteristic of granular filters and can be described by a CSD.Here, an extensive set of CSDs was calculated from data of 28 DEM simulations of varying coefficient of uniformity up to Cu = 6, PSD shape and relative density, using the numerical weighted Delaunay method proposed by Al-Raoush et al. (2003) and Reboul (2008).The user-defined merging overlap parameter used to join adjacent Delaunay tetrahedra was varied from 0% to 100%.As the merging overlap decreases, more tetrahedra merge to form larger voids and therefore fewer constrictions with a larger average diameter are identified.This highlights the fact that, for a granular material, void boundaries are subjective and no unique CSD can be identified.However, the trends observed for a single merging overlap value were the same regardless of whether the value was 0%, 50% or 100% and therefore general conclusions about the relationship between PSD, relative density and CSD can be drawn.The calculated CSDs were compared to two practical methods of filter design: (a) the inscribed circle method to analytically define a full CSD based on the largest circle which can be inscribed between combinations of three or four touching spheres (Locke et al., 2001) Fig. 9. Normalised CSDs obtained with weighted Delaunay method using a 50% overlap for loose, medium dense and dense and DEM samples with Cu = 1•5 and Cu = 6  (Kenney et al., 1985).Based on this analysis, the following conclusions can be drawn.
(a) For a given merging overlap value, the CSDs are similarly shaped regardless of the filter uniformity.
For linearly graded filters with the same smallest particle, D 0F , the CSDs become coarser as filter Cu increases from Cu = 1•2 to Cu = 3.For samples with Cu ! 3 the CSD curves form a narrow band.This shows very good qualitative agreement with the experimental work of Kenney et al. (1985), who found that, when normalised by D 0F , the controlling constriction size increased as filter Cu increased from Cu = 1•2 to Cu = 3, and then remained approximately constant up to at least Cu = 12.(b) Constriction sizes reduce with increasing relative density.However, the CSD shape remains similar.A difference in constriction diameter of around 10% was found between the CSDs from samples in the loosest and densest states analysed.
(c) A comparison between CSDs from the numerical weighted Delaunay and the analytical inscribed circle methods gave reasonable agreement for filters of low Cu values (Cu 2) and moderate relative densities in the inscribed circle method (R d = 30-60%).However, the large range of possible inscribed circle CSDs, depending on R d , was not replicated in the weighted Delaunay CSDs.Agreement between the methods became progressively worse as Cu increased to Cu = 6.
In particular the inscribed circle method resulted in much larger diameter constrictions.This poor agreement was due to an over-idealisation of soil fabric in the inscribed circle method.The assumption that the probability of a filter particle forming a constriction is proportional to surface area was found to be poor for more widely graded filters.This was due to coarse particles being separated from one another by smaller particles in these materials.Kenney et al. (1985) suggested were useful for practical purposes.In agreement with Kenney et al. (1985), D 5F was found to give narrowest band of CSD curves for filters with Cu ! 2, whereas the more commonly used D 15F was found to give the narrowest band for filters with Cu 3, regardless of whether the PSD was linear or bilinear.These findings give fundamental support to the use of characteristic diameters in filter design.
The CSDs presented here were calculated for homogeneous DEM samples of spherical particles.Future work should consider the fabric effects resulting from changes in particle shape, anisotropy and inhomogeneity (i.e.segregation) in real filters.In particular, micro-computed-tomography-based techniques (e.g.Taylor et al., 2015) allow the effect of particle morphology and soil fabric to be assessed at the microscale from three-dimensional images of real soils.

Fig. 3 .Fig. 4 .
Fig. 3. Two-dimensional analogy of three-dimensional merging of Delaunay cells to form voids (adapted from Al-Raoush et al., 2003): (a) correct constriction locations; (b) Delaunay triangulation of the void space; (c) false identification of constrictions based on Delaunay cell edges; (d) tangent sphere overlap used for Delaunay cell merging; (e) merging of Delaunay cells to form larger voids.Note that two-dimensional triangles represent three-dimensional tetrahedral; triangle sides represent tetrahedra faces and tangent circles represent tangent spheres

Fig. 5 .
Fig. 5. Comparison of CSDs calculated with the inscribed circle method and the weighted Delaunay (denoted WD) method: (a) Cu = 1•5; (b) Cu = 6.The three-and four-particle arrangements are shown in Figs 1(a) and 1(b) Fig. 6.(a) Cu = 1•5: comparison of diameters of filter particles making up vertices of constrictions according to inscribed circle and weighted Delaunay methods; (b) schematic diagram of soil fabric for uniform soils in which particles of all sizes can bound constrictions

Fig. 7 .
Fig. 7. (a) Cu = 6: comparison of diameters of filter particles making up vertices of constrictions according to inscribed circle and weighted Delaunay methods; (b) schematic diagram of soil fabric for widely graded soils in which larger particles are separated by smaller particles

Table 1 .
Summary of DEM simulations