Shape-dependent effective diffusivity in packings of hard cubes and cuboids compared with spheres and ellipsoids

We performed computational screening of effective diffusivity in different configurations of cubes and cuboids compared with spheres and ellipsoids. In total, more than 1500 structures are generated and screened for effective diffusivity. We studied simple cubic and face-centered cubic lattices of spheres and cubes, random configurations of cubes and spheres as a function of volume fraction and polydispersity, and finally random configurations of ellipsoids and cuboids as a function of shape. We found some interesting shape-dependent differences in behavior, elucidating the impact of shape on the targeted design of granular materials.


Introduction
Granular materials include a broad range of two-phase materials, both natural and synthetic, consisting of discrete solid particles (granules) surrounded by a continuous void phase. Some examples are columns for separation, chromatography, and catalytic reactions, cathode materials for lithium-ion batteries, and pharmaceutical materials for controlled release. Consequently, there is great interest in determining and understanding the effective macroscopic properties of these materials, such as effective diffusivity. [1][2][3][4][5][6] Understanding the impact of microstructural morphology is a key prerequisite for performing the targeted optimization of a material for a specific purpose. 7 The prototypical examples of granular materials are lattices and random configurations of impermeable spheres. However, contemporary particle synthesis techniques facilitate the production of a broad range of anisotropic nanoparticles, shaped like ellipsoids, rods, disks, cubes, icosahedra, tetrahedra, triangles, prisms, half shells, tripods, and stars, to name a few. [8][9][10][11][12][13][14][15][16][17][18][19][20][21][22] The realization that particles with anisotropic shapes (and also anisotropic interactions) can be used as colloidal building blocks to yield complex structures with novel properties has attracted substantial attention. [23][24][25][26] A vast body of literature exists on the topic of understanding shape-dependent assemblies of particles including their phase behavior, phase transitions, random close packing, and spontaneous shape-induced crystallization. Just to give a brief account, studies performed include the random close packing of ellipsoids, 27 monodisperse and polydisperse hard spheres, 28,29 polyhedra including platonic solids, [30][31][32][33] superballs and superellipsoids, [34][35][36] cubes with round edges, 37 and monodisperse and polydisperse cubes and cuboids. [38][39][40][41][42] Hard cubes and cuboids and their packing and phase behavior have been studied both as limiting cases of superballs and superellipsoids, 35,36 and modeled exactly as cubes and cuboids. [39][40][41] The true nature of their behavior is somewhat elusive; as stated by Jiao and Torquato, 31 attempts to create random close packings (or rather maximally random jammed packings) of cubes easily lead to high degrees of order, raising questions concerning the appropriateness of some of the algorithms suggested so far to generate these packings. Nonetheless, random sequential addition [39][40][41] can produce packings with a high degree of randomness, as can the order-constrained stochastic optimization method described in ref. 42.
Effective transport coefficients in random, two-phase materials, and in particular the case of spherical granules is a rather well investigated subject including both effective diffusivity (including other physical processes mathematically analogous to effective diffusivity, such as electrical conductivity, thermal conductivity, and magnetic permeability) [43][44][45][46][47][48][49][50][51][52] and fluid permeability, [53][54][55][56][57][58] oftentimes incorporating microstructural descriptors such as the n-point correlation (probability) functions for n = 1, 2,. . . introduced by Brown. 59 In light of the increasing capabilities to manufacture anisotropic nanoparticles, the topic of effective macroscopic properties and the impact of shape thereon has become increasingly interesting. Transport properties with anisotropic granules have been investigated in several studies. In ref. 60, random packings of granules shaped as spheres, ellipsoids, cylinders and parallelepipeds were generated using sequential deposition in a gravitational field, and their conductivity, permeability, and diffusivity were studied (but the latter only using ellipsoids for computational reasons). In ref. 61, the impact of shape on the effective permittivity of ellipsoids was studied. In ref. 62, different irregular particle shapes constructed as the union of overlapping spheres were created and the impact on permeability was investigated. Nonetheless, not so much attention has been directed toward aspherical particles and effective diffusivity.
In this work, we performed computational screening of effective diffusivity in packings of (hard and impermeable) cubes and cuboids. Very little work is done in terms of both analytical models and simulations of effective diffusivity for particles with these shapes. As a reference, we compared all the results with packings of spheres and ellipsoids. We studied both periodic (lattice) configurations, i.e. simple cubic crystals and face-centered cubic crystals, and random configurations for different solid volume fractions, polydispersities and shapes. In total, more than 1500 structures are generated and screened for effective diffusivity. The aim of this computational screening paradigm is to explore the effect of varying different geometrical parameters independently, and discover and understand generic design rules for the tailoring of effective transport properties. This effort will aid in guiding future experimental work.

Structure generation
We first describe the details of generating the lattice structures in a cubic, periodic simulation domain with side L. A simple cubic (SC) lattice of spheres is generated by placing a single sphere with radius r = 1 in (L/2,L/2,L/2) with L being defined by the solid volume fraction f through The upper bound for the solid volume fraction is f = p/6 E 0.5236, which is when the spheres tangent each other. An SC lattice of cubes is generated by placing a single cuboid with semi-axes r 1 = r 2 = r 3 = r = 1 in the same position, with L being The upper bound for the solid volume fraction is f = 1 when the cubes tangent each other. A face-centered cubic (FCC) lattice of spheres is generated by placing four spheres with radii r = 1 in (0,0,0), (L/2,L/2,0), (L/2,0,L/2), and (0,L/2,L/2), with L being The upper bound for the solid volume fraction is f ¼ p= 3 ffiffi ffi 2 p À Á % 0:7405 when the spheres tangent each other. An FCC lattice of cubes is generated by placing four cuboids in the same positions as above, with L being The upper bound for the solid volume fraction is f = 1/2 when the cube corners tangent each other, resulting in a 3D 'checkerboard' structure. The random microstructures are generated using in-house developed software implemented in Julia (www.julialang.org) 63 and available in a Github repository (https://github.com/rod ing/whitefish_generation, version 0.2). Random configurations of non-overlapping particles are generated using a hard particle Markov Chain Monte Carlo (MCMC) algorithm. For spheres, the overlap condition is straightforward. For ellipsoids, the Perram-Wertheim criterion 64 for two ellipsoids of arbitrary orientation is used. For cubes and cuboids, a separating axis theorem 65 is employed for the detection of overlap for two cuboids. To quantify the degree of overlap, it is further determined how far one of the particles needs to move along the axis defined by their respective center points to remove the overlap. First, particles are assigned uniformly distributed locations and orientations (the latter encoded using a quaternion representation). Second, the configurations are relaxed by sequentially performing random translations of all particles and then random rotations of all particles until no two particles overlap. Proposed translations and rotations are only accepted if they lead to a lower or equal degree of overlap for the considered particle. These 'local' stochastic optimization steps eventually lead to a 'global' optimization, resulting in no overlap. Third, the configurations are equilibrated by performing a large number of random translations and rotations ensuring a distribution in location and orientation that is as uniform as possible. Now, if the desired solid volume fraction f is less than 0.50, we are done. Otherwise, as a final step, the steps above are performed for f = 0.50 and the configuration is compressed in small steps, Df = 10 À5 , until the target solid volume fraction f target is reached. The proposed translations are normal distributed with standard deviation s t in each direction. The proposed rotations are normally distributed with standard deviation s r in a random direction. In every step, s t and s r are chosen in an adaptive fashion to aim for an acceptance probability of 0.25. For very high solid volume fractions, the configuration might get jammed before reaching f target . In this case, the generation is restarted from scratch with a new random configuration. All configurations have 512 particles (there is a point to having an even cube, i.e. 8 3 , because otherwise the configurations of the cubes actually cannot reach high solid volume fractions in the (artificially) periodic simulation domain). Using a dual Intel Xeon E5-2699 v4 setup, the execution times for the generation are between B1 min and B2.5 hours (single thread), depending strongly of course on the type of particle and on f target .

Diffusion simulation
Obstructed diffusion in the generated structures is simulated using an in-house developed software also implemented in Julia (www.julialang.org) 63 in a parallel implementation and also available in a Github repository (https://github.com/rod ing/whitefish_diffusion, version 0.2). A 'random walk particle tracking' technique 50 is used, initially placing an ensemble of N = 10 6 point particles uniformly distributed in the void part of the simulation domain, by rejection sampling. The stochastic particle motion is simulated as a Gaussian random walk with time resolution dt = 2.5 Â 10 À5 (a.u.) and corresponds to the diffusion coefficient D 0 = 1 (a.u.). Hence, random normal distributed displacements are added to the current position in each time step with zero mean and standard deviation ffiffiffiffiffiffiffiffiffiffiffiffi ffi 2D 0 dt p . The positions are recorded at each major time step Dt = 0.01 (a.u.). The simulation proceeds up to t max = 50 (a.u.) or 100 (a.u.) (for different data sets). Single and multiple rejection boundary conditions are implemented but multiple rejection is used here; 66 proposal displacements are generated until one is found that does not mean leaving the void phase. A cell list data structure (with 10 Â 10 Â 10 cells for large random configurations, but only a single cell for lattice structures) is used to accelerate computations. The computations are run on a variety of architectures comprising around 20 cores yielding execution times between B2 hours and B30 hours, within the total 5 Â 10 5 core hours used.
A time-dependent effective diffusion coefficient (i.e. obstruction factor) is computed as x n ðtÞ À x n ð0Þ ð Þ 2 þ y n ðtÞ À y n ð0Þ ð Þ 2 þ z n ðtÞ À z n ð0Þ ð Þ 2 h i : Here, x n (t), y n (t), and z n (t) are the positions of the n:th particle at time t. The effective diffusivity is obtained as the asymptotic value D N /D 0 A (0,1). In Fig. 1 some examples of computed effective diffusivity curves are shown.

Results and discussion
First, we consider the lattice structures, i.e. simple cubic (SC) and face-centered cubic (FCC) configurations of spheres and cubes. Examples of the structures and the results of the diffusion simulations are shown in Fig. 2. For the spheres, the results are in very good agreement with earlier analytical findings 67 (see Fig. 2 in this paper, and the solid black lines in the figure herein) and with earlier simulated findings using a very similar simulation technique as herein. 50 Obviously, for small f the impact of shape is small, and grows larger for increasing f. The FCC of cubes consistently has the lowest diffusivity of them all. For f smaller than B0.4, the SC of cubes has smaller diffusivity than the SC of spheres, whereas for f larger than B0.4, the opposite is true. Similarly, for f smaller than B0.6, the SC of cubes has a smaller diffusivity than the FCC of spheres, whereas for f larger  than B0.6, the opposite is true. The 3D checkerboard structure of the FCC of cubes yields very narrow 'channels' between large open spaces, even for relatively low solid volume fractions; indeed, for f = 0.5, the void phase is discontinuous and the effective diffusivity is zero. Furthermore, there is some benefit in terms of effective diffusivity in the SC and FCC of spheres as compared to the FCC of cubes as long as the solid volume fraction is not too large, because of a less tortuous path for the diffusing species with smoothed boundaries of the void phase. However, for large solid volume fractions, the SC of cubes can provide highly non-tortuous and straight paths (at least along the axes), simply because packing of cubes is more efficient than that of spheres, leading to rather high diffusivities even for extremely high solid volume fractions. We first believed this to be an artifact of our simulation method, in particular the multiple-rejection boundary conditions, but as can be seen below in Fig. 3, when the 'straight channels' are no longer there, the diffusivity indeed drops to zero much more quickly (the curve shown in Fig. 2 Fig. 3. For low solid volume fractions, spheres yield a higher effective diffusivity than cubes, the same as above. For high solid volume fractions, it appears that the benefit of straight 'channels' induced by the cubic shape in lattices is lost for random configurations. This is not surprising; rather, it is intuitively obvious that adjacent cubes randomly displaced from each other can provide a very efficient obstacle. For f larger than B0.6, we observe akin to ref. [39][40][41][42] that ordering and orientational correlations start manifesting themselves (for f Z 0.57, crystallization behaviour has been observed). It can hence be discussed whether these configurations are random and to what extent; herein, it suffices to conclude that our algorithm will create crystals for high solid volume fractions. The substantial spread in effective diffusivity for cubes with solid volume fractions larger than B0.6 indicates that it is largely up to chance whether the configuration will provide 'simple' paths for the diffusers or not. This may also be highly dependent on the choice of algorithm for generating the configurations. We have a hypothesis that the slightly slowed decrease (the slight 'bump') in diffusivity in the range f = 0.6-0.8 depends on the increased ordering among the cubes yielding some 'channels' between the cubes, before the effect of a larger solid volume fraction once again starts to dominate. In any case, a more detailed answer on the true effective diffusivity for cubes with very large solid volume fractions would likely require much larger configurations.
Third, we investigate the effect of polydispersity for f = 0.5 by introducing a lognormal distribution for the particle volumes (due to the family of lognormal distributions being closed under powers, the semi-axes/radii are also lognormally distributed). Examples of the structures and the results of the diffusion simulations are shown in Fig. 4. The coefficient of variation (CV), i.e. the ratio of standard deviation to the mean, is varied from 0 to 3. The effect of polydispersity is not that pronounced, but a slightly decreasing trend for increasing the coefficient of variation is seen. As the calculation of effective diffusivity is a scale-independent problem, increasing the size of all particles would not change anything. However, introducing the occasional large particle through polydispersity will apparently contribute toward blocking the diffusion more efficiently. Although we only investigate for f = 0.5, it is clear that the effect of polydispersity will be positively correlated with solid volume fraction and vanish for f approaching zero. Fourth, we investigated the effect of shape for ellipsoids and cuboids for f = 0.5. For both particle types, we studied semi-axes 1 : 1 : a, where the semi-axis ratio a is varied between 0.25 and 4. Examples of the structures and the results of diffusion simulations are shown in Fig. 5. As can be expected from the results for both mono-and polydisperse packings of spheres and cubes, cuboids do consistently yield lower effective diffusivities than ellipsoids. However, for very oblate particles (to the extent that this is an appropriate word for cuboids), the difference gradually vanishes. The same convergence is not observed for very prolate particles in the investigated a regime, although it seems plausible that for extremely elongated particles, shape (of the cross-section) would once again not matter so much. We investigated this for a = 10 and f = 0.4 (because we cannot reach f = 0.5 for such elongated particles) and somewhat surprisingly found that D/D 0 was B0.77 for ellipsoids and B0.73 for cuboids, i.e. the difference in diffusivity reduces from B0.05 to B0.04 when jumping from a = 4 to a = 10, suggesting that they 'converge' very slowly as a increases. We cannot generate structures with arbitrary large values of a due to numerical issues with the overlap criterion (for ellipsoids) and particles possibly overlapping themselves due to the periodicity of the simulation domain, so we have to settle for this result. Anyhow, the two particle types have the same generic relation between shape and effective diffusivity, suggesting some degree of universality with respect to the aspect ratio impact on diffusivity. The importance of shape and aspect ratio could be due to the shape itself, and also due to the ordering that the shape can impose on the overall structure, or a combination of the two. A deeper analysis of microstructural descriptors would be required to find an answer. Although we only investigated for f = 0.5, it is clear that the effect of shape will be  positively correlated with solid volume fraction and vanish for f approaching zero.

Conclusion
We have performed computational screening of effective diffusivity in lattices and random configurations of cubes and cuboids compared with spheres and ellipsoids. In total, more than 1500 structures were generated and screened for effective diffusivity. Very little is done in terms of both analytical models and simulations of effective diffusivity in the literature for cubic and cuboidal particles. For lattice configurations, the solid volume fraction will determine which shape gives higher and lower effective diffusivity. For random configurations, cubes and cuboids consistently yield lower effective diffusivities. Although the solid volume fraction is the primary determinant of effective diffusivity, shape effects are still substantial and more research is needed to understand them. The shapedependent differences found in the behavior help in understanding the impact of shape on the targeted design of granular materials, with respect to the relations between porosity, specific surface, pore size distributions, and effective properties. Further work would naturally include an extension to superballs/superellipsoids as well as cubes/cuboids with rounded edges and corners, both of which can 'interpolate' between spheres/ellipsoids on the one hand and cubes/cuboids on the other. Furthermore, the use of 2-and 3-point correlation functions as well as other microstructural descriptors will help elucidate the nature of the shape dependence in these systems.

Conflicts of interest
There are no conflicts to declare.