Generation of Granular Media

Transport in Porous Media 26: 99–107, 1997. c 1997 Kluwer Academic Publishers. Printed in the Netherlands. 99 Generation of Granular Media LAURENT TACHER, PIERRE PERROCHET and AURÈLE PARRIAUX Laboratory of Geology, Swiss Institute of Technology, GEOLEP-EPFL, CH-1015 Lausanne, Switzerland; e-mail: [email protected] (Received: 5 March 1996; in final form: 24 September 1996) Abstract. A discrete reduced distance method to generate 2-D and 3-D granular porous media is presented. The main property of the method is to produce heterogeneous and/or anisotropic packed beds of joined grains with arbitrary shapes and optimum fitting (i.e., minimum porosity). The iterative generation process starts with the coarsest grain and adjusts the size and location of the next ones depending on the updated available space. Hence, grain size distribution cannot be specified directly but is merely the consequence of user defined input parameters. The latter consists of a set of randomly distributed initial points, a few typical predefined grain shapes as well as the minimum and maximum grain diameters. The simulated granular media can readily be processed by an appropriate mesh generator to allow for subsequent numerical solutions of differential equations. Key words: modelling, microscopic scale, granular media. 1. Introduction The derivation of classical macroscopic flow and transport parameters – such as permeability, solute and thermal dispersivities, sorption isotherms, reaction and decay coefficients, etc. – based on an integration of pore scale processes is a more than ever open problem (for review article see Dullien, 1991; Adler and Thovert, 1993; Allaire, 1989; Barrere et al. 1992; Blunt and King, 1990; Borne, 1992; Brenner, 1980; Chapman, 1992; Palaniappan, 1993; Withaker, 1986. In numerous and recent studies, promising attempts are made to relate permeability and dispersion to their microscopic origins. In most cases this is done by solving local or average field equations at the pore level of reconstructed porous structures ranging from the basic orthogonal assembly of homogenous spheres (e.g., German, 1989; Eidsath et al., 1983) to RUCs (Prieur du Plessis and Masliyah, 1991), interconnected networks (Blunt and Kin, 1990), fractal carpets (Adler, 1988; Adler and Jacquin, 1987) and random media with autocorrelated porosity (Sallès et al., 1993). These artificial geometrical structures aim at simulating relevant microscopic features (fluid/solid area of contact, tortuosity) governing the values of the macroscopic parameters. Other studies consider the behaviour of packing under the action of various physical forces (Jullien and Meakin, 1992; Yen and Chaki, 1992; Békri et al., 1995), but starting from a structured or random grain distribution. The grains’ shapes are in most cases circular/spherical or elliptical (Thies-Weesie et al., 1995). 100 LAURENT TACHER ET AL. Figure 1. Joined grains generation stages: (a) set of random IP s, (b) first grain centre C1 , (c) final granular bed. However, due to the difficulties inherent in generating polydispersed granulometries and meshing the grains or the voids between them in two and three dimensions, computerization of realistic granular media including individual grain details does not seem to have been performed yet. This paper presents the basic principles of a relatively simple method developed to take into account these microscopic geometrical details in 2- and 3-D. Unlike the approaches using fractals or autocorrelation functions, the shapes of the grains (possibly random) are predefined. The generation process is described and illustrated by some typical examples. 2. Principles In the domain considered, a starting set of random initial points (IP s) is produced using a density function representative of the largest grains to be created (Figure 1(a)). A first grain, with centre C1 , is then positioned at that location which maximizes the shortest Euclidean distance between C1 and any IP s. In other words C1 is positioned where it is as ‘lonely’ as possible, and the distance to its closest neighbour is the radius of the grain (Figure 1(b)). The points belonging to that grain surface are then added to the IP list and a second grain centre C2 is looked for in the same manner, but with the updated IP list. This simple procedure is carried out (Figure 1(c)) until a specified smallest grain radius is reached. For example, if the minimum and maximum radii are equal, a monodisperse packing is built. At any stage of the generation, one notes that each grain (circles in 2-D or spheres in 3-D) satisfies the Delaunay criterion (Weatherill, 1992), which means that it is bounded by at least three or four of the IP s in 2-D or 3-D respectively, and that there are no IP s inside it. The only parameters to provide are the mean distance between the IP s before the first iteration (which is also the maximum grain diameter), and the minimum grain radius, a porosity or a total number of grains. Both criteria may be used to stop the process. Since the first few grains match exactly some of the input random initial points, their radius may slightly exceed the required maximum radius. Thus, 101 GENERATION OF GRANULAR MEDIA to respect strictly the maximum radius, the radii of the grains can be reduced as follows: Cn = Minimum(Cn ; Cmax ): Practically, the search for every Cn locations is performed by means of a regular orthogonal background grid G, the size of which is finer than the smallest grain. A large number of possible grain centres are thus produced (G) and then, at each iteration n, the specific location Cn showing the maximum distance d to the closest I P is sorted out following the rule Cn 2G such that Minimum(d(Cn ; I P(1;jn ) )) = Maximum over G; (1) where jn is the number of initial points at the start of iteration n (number of random initial input points plus points belonging to the surface of the n 1 previously generated grains). After Cn is selected, it is necessary to update the distance of all points g of to the initial points, now including the new points belonging to the surface of the grain centred at Cn . Since these distances need to be modified only by the last generated grain, this update is the single minimum operation G (n) d (g; I P(1;jn ) ) = Minimum(d(n 1) (g; I P(1;jn 1) (n) ); d (g; I P(jn 1 ;jn ) )); (2) where I P(jn 1 ;jn ) are the points of the grain found at iteration n, and added in the I P list. Consequently the above process is quite fast and becomes particularly efficient in 3-D. As an example the porous medium illustrated on Figure 2, involving a 100  100  100 grid G with minimum and maximum radii of 1.5 and 20 respectively, is built in five minutes with a 100 MHz processor. The total number of grains is 1714 and the porosity is 25.8 per cent. At each iteration the size of the new grain is maximum with respect to the available space, so that the process fills the domain with gradually decreasing size grains and optimum fitting. Thus, the final grain size distribution is not directly imposed but is a consequence of the given parameters. However, one can limit not only the maximum radius, as above, but all created grains, according to a given distribution, which leads to leaving room between them. One has the choice between:  not limiting the radii according to a distribution and thus getting a realistic stacking, since each grain is connected to at least three of its neighbours,  limiting the radii to respect a distribution, which may leave room between the grains. 102 LAURENT TACHER ET AL. Figure 2. 3-D porous medium model with 1714 spherical grains. Porosity = 25 8 per cent. : At this stage, the distances are Euclidean distances so that only circular or spherical grains can be stacked. To address elliptical and other grain shape fitting, it is obviously possible to create a set of circular/spherical grains and to replace each of them with a contained shape, but this does not lead to a realistic fitting. Thus, it is necessary to introduce reduced distances and discrete reduced distances. This means that the distances are to be measured differently according to the direction. 3. Elliptic Grains Distances of points g of G to an initial point IPj are reduced distances d , illustrated on Figure 3 and defined in 2-D according to d (g; IP  j )= " 1 0 0 1=a #" cos cos sin cos #(
x
y #) ; (3) where a is the flatness coefficient. p Substitution of d x2 yields: =
+
y2 in (1) by d from the above equation simply C n 2G (d (C such that Minimum n ; IP(1 ;j n) )) = Maximum over G; (4) GENERATION OF GRANULAR MEDIA 103 Figure 3. Sketch to calculate the reduced distance d for elliptical grains. Figure 4. 2-D porous medium model with elliptical grains and horizontal anisotropy. and the update of the distances of the points g of G to initial points is performed using the same reduced distances. In addition, grain flatness (stretching) and orientation may vary randomly at each iteration to simulate grains of heterogeneous shapes, sizes and orientations. When a given anisotropy is desired, the angle is set accordingly. The value = 0 corresponds to a horizontal stratification. This is illustrated on Figure 4 where 140 elliptical grains are generated on a 400  400 grid G. The minimum and maximum 104 LAURENT TACHER ET AL. Figure 5. Sketch to calculate the discrete reduced distances d i for grains of arbitrary shapes. horizontal radii are 7 and 70 respectively and the flatness is at random between 1.5 and 6. 4. Arbitrary Shapes More generally, arbitrary shapes may be created and mixed by using discrete reduced distances. Each grain is defined by a set of discrete reduced distances  di , which are the distances of its surface points to Cn in a reference local space (Figure 5). For a given grain, the location and size providing the optimum fitting is searched as previously by Cn 2G such that Minimum(d (Cn ; I P(1;jn ) )) = Maximum over G; i (5) where subscript i varies from 1 to the number of points required to describe the shape of the grain (e.g., 35 points). At each iteration, the shape of the grain may be chosen randomly among a predefined set. Furthermore, each predefined shape can be modified at random. An example is given on Figure 6 where the grid G is 150  150 and the size of the generated 209 grains lies between 1 and 30. In this particular case, two types of shapes are selected in equal proportions but with random orientations. One of the types is a circle defined with 35 points, with d lying randomly between 0.9 and i 1 in a reference space. The other one is constant. 5. Pore Scale Analyses With the kind of results presented above at hand, a few issues can directly be addressed without further difficulties. Actually, the generation of realistic granular media together with all relevant microscopic details make it straightforward to analyse possible statistical relationships between, say, porosity or specific contact area and grain size (and shape) distributions. It is indeed well known that these two key-parameters largely determine the conductive properties of granular packed beds. GENERATION OF GRANULAR MEDIA 105 Figure 6. 2-D porous medium model with arbitrary grains and random rotations. As far as dynamic processes such as flow, transport or deformation are concerned, it is clear that the domain discretization step must be performed before any numerical scheme is used. Depending on what is to be simulated, one may need to discretize either the solid or the porous phase, or both. Here, the mesh generation technique of Tacher and Parriaux (1996) can be used to fill the pore volume with triangular elements (tetrahedra are used in 3-D). This technique was specially designed to handle domains of complex geometries and works basically with exactly the same principles as the grain generator. Figure 7 illustrates the resulting computational mesh consisting of 3557 nodes and 4950 elements. The reconstructed granular medium is the same as in Figure 4 but the smallest diameters are filtered out to preserve a good visualisation of the triangles. 6. Final Remarks The method presented here allows for a flexible generation and a detailed description of virtually any granular packed bed. It seems to be realistic enough to simulate a wide range of sedimentary media where porosities are reasonably minimized. However, for magmatic and metamorphic rocks, the choice of grain shapes at each iteration should rather be made with respect to the remaining available space than predefined, in order to reproduce fitting grains or crystal-like patterns. 106 LAURENT TACHER ET AL. Figure 7. Triangular mesh of a reconstructed granular medium pore volume. In association with an appropriate meshing technique the present method can be a tool to start a series of new numerical experiments regarding the derivation of macroscopic parameters. Actually, with accurate meshes of both the grains and the pore volume, one can expect to gain further relevant insights in flow, thermal, chemical or mechanical issues using deterministic modelling approaches. This can be done either by solving the primitive variables of local equations (e.g., Stokes equations, advection equations, diffusion equations) or by tackling various closure schemes as stated by Brenner (1980), Whitaker (1986), Barrère et al. (1992) regarding dispersion and permeability tensors. Other applications may deal for instance with the effects of the grain stacking on porosity or with the sensitivity of the porosity/granulometry relationship to the effective granular structure of the medium. GENERATION OF GRANULAR MEDIA 107 References Adler, P. M. and Thovert, J.-F.: 1993, Fractal porous media, Transport in Porous Media 13, 41–78. Adler, P. M. and Jacquin, C. G.: 1987, Fractal porous media I: Longitudinal Stokes flow in random carpets, Transport in Porous Media 2, 553–569. Adler, P. M.: 1988, Fractal porous media III: Transversal Stokes flow through random and Sierpinski carpets, Transport in Porous Media 3, 185–198. Allaire, G.: 1989, Homogenization of the Stokes flow in a connected porous medium, Asympt. Anal. 2, 203–222. Barrere, J., Gipouloux, O. and Whitaker, S.: 1992, On the closure problem for Darcy’s law, Transport in Porous Media 7, 209–222. Békri, J., Thovert, J. F. and Adler, P. M., 1995, Dissolution of porous media, Chem. Engng. Sci. 50, 2765–2791. Blunt, M. and King, P.: 1990, Macroscopic parameters from simulations of pore scale flow, Phys. Rev. 42(8), 4780–4787. Borne, L.: 1992, Harmonic stokes flow through periodic porous media: a 3D boundary element method, J. Comput. Phys. 99, 214–232. Brenner, H.: 1980, Dispersion resulting from flow through spatially periodic porous media, Philos. Trans. R. Soc., London, Ser. A 297, 81–133. Chapman, A. M. and Higdon, J. J. L.: 1992, Oscillatory Stokes flow in periodic porous media, Phys. Fluids A, 4, 2099–2116. Dullien, F. A. L.: 1991, Characterisation of porous media – Pore level, Transport in Porous Media 6, 581–606. Eidsath, R., Carbonell, R. G., Withaker, S. and Hermann, L. R.: 1983, Dispersion in pulsed systems III: Comparison between theory and experiments for packed beds, Chem. Engng. Sci. 38, 1803–1816. German, R. M.: 1989, Particle Packing Characteristics, Metal Powder Industries Federation, New York, 433 pp. Jacquin, C. G. and Adler, P. M.: 1987, Fractal porous media II: Geometry of porous geological structures, Transport in Porous Media 2, 571–596. Jullien, R., Meakin, P. and Pavlovitch, A.: 1992, Three-dimensional model for particle-size segregation by shaking, Phys. Rev. Lett. 69(4), 640–643. Palaniappan, D., 1992, Arbitrary Stokes flow past a porous sphere, Mech. Res. Comm. 20, 309–317. Prieur du Plessis, J. and Masliyah, J. H.: 1988, Mathematical modelling of flow through consolidated isotropic porous media, Transport in Porous Media 3, 145–161. Prieur du Plessis, J. and Masliyah, J. H.: 1991, Flow throught isotropic granular porous media, Transport in Porous Media 6, 207–221. Sallès, J., Thovert, J. F. and Adler, P. M.: 1993, Reconstructed porous media and their application to fluid flow and solute transport, J. Contam. Hydrol. 13, 3–22. Tacher, L. and Parriaux, A.: 1996, Automatic nodes generation in a N-dimensional space. Comm. Num. Meth. Engng. 12(4), 243–248. Thies-Weesie, D. M. E., Philipse, A. P. and Kluijtmans, S. G. J. M.: 1995, Preparation of sterically stabilized silica-hematite ellipsoids: Sedimentation, permeation, and packing properties of prolate colloids, J. Colloid Interface Sci. 174, 211–223. Withaker, S.: 1986, Flow in porous media I: A theoretical derivation of Darcy’s law, Transport in Porous Media 1, 3–25. Weatherill, N. P.: 1992, The Delaunay triangulation in CFD, Comput. Math. Appl. 24(5/6), 129–150. Yen, K. Z. Y. and Chaki, T. K.: 1992, A dynamic simulation of particle rearrangement in powder packings with realistic interactions. J. Appl. Phys., 71–7, 3164–3173.