Eﬀective diﬀusivity in lattices of impermeable superballs

,

For exact prediction of effective properties in multi-phase media, including diffusivity, complete microstructural information is required.However, analytical expressions for diffusivity are scarce even when the microstructure is known and can be described in an analytic fashion.One exception is that of lattices of permeable and impermeable spheres (impermeable implying that the diffusion takes place only in the intermediate space outside the spheres), where analytical expansions have been derived [33].A useful tool for the characterization of microstructures is the n-point probability/correlation functions introduced by Brown [34].The correlation functions for all n combined correspond to complete knowledge of the microstructure; however, in reality they can be computed only for n ≤ 3. Consequently, deriving bounds and approximations based on lower-order correlation functions has attracted substantial interest.Optimal two-point bounds i.e. involving correlation functions for n ≤ 2 were derived by Hashin and Shtrikman for effective magnetic permeability in isotropic two-phase media [35], a case which is mathematically analogous to effective diffusivity.These two-point bounds were improved upon by the three-point bounds introduced by Beran [36], which were in turn simplified by Torquato [37] and Milton [38] to involve the solid volume fraction and a three-point parameter known as ζ 2 .
In previous work [39], we investigated the shape-dependent effective diffusivity in packings of cubes and cuboids compared with spheres and ellipsoids.This study elucidated some interesting aspects of the impact of shape in both random configurations and lattices.Elaborating on that work, we herein study lattices of impermeable superballs, which are an 'interpolating' class of shapes that include both spheres and cubes as special cases.Superballs (and their two-dimensional counterpart, superdisks) were first introduced and studied in the context of granular media by Torquato and colleagues [40,41].We perform stochastic simulations of diffusion, where point particles diffuse in the intermediate space between the solid superballs, to estimate effective diffusivity as a function of shape and solid volume fraction for simple cubic (SC), body-centered cubic (BCC), and face-centered cubic (FCC) lattices, and study the performance of an analytical approximation based on twoand three-point correlation functions and the microstructural parameter ζ 2 .

II. DIFFUSION SIMULATION
Obstructed diffusion is simulated using in-house software implemented in Julia (www. julialang.org)[42], available in a Github repository (https://github.com/roding/whitefish_diffusion, version 0.3).Random walk-based simulation [43,44] is used.Initially, an ensemble of N = 4 • 10 6 point particles (i.e.infinitesimally small particles) are distributed uniformly in the void phase of a periodic, cubic simulation domain [0, L] 3 (the value of L depends on the type of lattice as well as on the desired solid volume fraction, and will be elaborated upon below).Particle motion is simulated as a Gaussian random walk with time resolution δt = 2.5 • 10 −5 (a.u.) and corresponding to diffusion coefficient D 0 = 1 (a.u.):Hence, random normal distributed displacements with zero mean and standard deviation √ 2D 0 δt are added to the current position in each time step.The positions are recorded at each major time step ∆t = 0.01 (a.u.).The simulation proceeds up to t max = 100 (a.u.).
Multiple rejection boundary conditions [45] are used, meaning that proposal displacements are generated until a feasible jump is found i.e. a jump such that the particle still resides in the void phase and does not move into any of the impermeable superballs.
A time-dependent effective diffusion coefficient (i.e.obstruction factor) is computed as 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 ∞ /D 0 ∈ (0, 1).In Fig. 1, some examples of computed effective diffusivity curves are shown.All the diffusion simulations use

III. CHARACTERIZATION AND PREDICTION
The two-point bounds for effective diffusivity in a solid-pore structure are where ϕ is the solid volume fraction, which is not very restrictive.The three-point bounds as derived by Torquato [37] and Milton [38] are where ζ 2 is microstructural parameter defined as Here, S n are the n-point correlation functions for n = 1, 2, 3 and are defined in the following way: S 1 is the probability that a random point is in the void space, S 2 (r) is the probability that two random points separated by a distance r are both in the void space, and S 3 (r, s, θ) is the probability that three random points forming a triangle with two sides r and s and with an angle θ between them are all in the void space.The one-point correlation function is trivially equal to 1 − ϕ.The higher-order functions are generally not as simple.The microstructural parameter ζ 2 can also be used for a three-point approximate analytical expression for the effective diffusivity, originally derived for the case of electrical conductivity by Torquato [46] and used in the modified form above for the case of diffusivity in several papers [43,44,47].Both the three-point bound and the three-point approximation approach the two-point bound as ζ 2 approaches zero.
We wish to employ Eq. ( 5) for prediction of effective diffusivity and to compare with simulation results.Earlier calculations of ζ 2 have typically been performed using a Monte Carlo integration scheme, where numerical integration is carried out over a circular grid, equidistant in θ and log-equidistant in r and s.The integral has been evaluated as the average of a large number of such numerical integrations using random positions and orientations for the grids [44].However, we propose to perform Monte Carlo integration using random sampling in the whole (continuous) domain r min ≤ r, s, ≤ r max and 0 ≤ θ ≤ π, approximating the integral as an average of N MC values of the integrand.The angle θ is sampled from a uniform distribution.To avoid the diverging factors 1/r and 1/s, we use importance sampling, generating values of r and s distributed according to a probability distribution f (r) proportional to these factors i.e.
Sampling from this distribution is performed by sampling log(r) and log(s) uniformly from Finally, ζ 2 is computed as Throughout this study, we let N MC = 10 12 .Also, we let r min = 10 −5 and r max = 10 2 .It was found that this choice captures the asymptotic properties of the correlation functions (data not shown), namely that and the last with t being the third side of the triangle, i.e. requiring that θ > 0 (for θ = 0 the limit as r and s approach ∞ would instead be (1 − ϕ) 2 ).The computation is implemented in Julia (www.julialang.org)[42] and is validated on sphere lattices, where the obtained values of ζ 2 are in very good agreement with those found before [43,44].All the computations use approximately 10 5 core hours.

IV. RESULTS AND DISCUSSION
A superball is a generalization of a sphere, defined by the equation where (x 0 , y 0 , z 0 ) are the center coordinates and r is the radius/semi-axis.For p = 2, a sphere with radius r is obtained and for p = ∞, a cube with sides 2r is obtained.We investigate the range 2 ≤ p ≤ ∞ here, which includes superballs with cube-like, rounded shapes (for p < 2, the shape is octahedron-like instead, but we do not investigate this regime).The volume of a superball is where β(x, y) = Γ(x)Γ(y)/Γ(x, y) is the beta function.This expression collapses to the expressions for sphere volume and cube volume in the limits of the range.We let r = 1 in all cases.
We investigate effective diffusivity for SC, BCC, and FCC lattices, and we do so for a range of p values for constant solid volume fractions ϕ.For SC lattices, we use ϕ = 0.40, 0.45, and π/6, the latter being the largest value attainable for all p, for p = 2, 3, 4, 5, 10, 20, and ∞.In Fig. 4, the simulated effective diffusivity (taken as the asymptotic value D ∞ /D 0 ) and the estimated effective diffusivity using Eq. ( 5) are shown.We also validate by comparing to the analytical results of Blees and Leyte [33] for the case of spheres, finding good agreement (a more detailed comparison of our simulation framework with their analytical model is performed in [39]).Generally, the estimate is better for lower solid volume fractions (higher effective diffusivity).Also, the estimate is better for p roughly between 3 and 10; for p = 2 and p > 10 the discrepancy increases.However, the estimates suggest that the ζ 2 parameter captures the impact of shape rather accurately.Additionally, we notice that there is an optimal value of p somewhere around 5 in terms of effective diffusivity, which we will further investigate below.For BCC lattices, we perform the same investigation, but for ϕ = 0.20 and FCC (red; dark gray).For p = 2, the well-known values ∼ 0.5236 (SC), ∼ 0.6802 (BCC), and ∼ 0.7405 (FCC) are attained; as p approaches infinity, the values 1(SC), 1/2 (BCC), and 1/4 (FCC) are attained.Notably, ϕ max is an increasing function of p for SC, but decreasing for BCC and FCC. and 0.25, the latter being the largest value attainable for all p.In Fig. 5, the simulated effective diffusivity and the estimate are shown.Also here, we find good agreement with the analytical results of Blees and Leyte [33] for the case of spheres.Interestingly, the estimate is better for smaller p, and somewhat better for smaller ϕ.Also, the effective diffusivity is monotonically decreasing as a function of p, in contrast to the case of SC lattices.Finally, for FCC lattices, the investigation is performed for ϕ = 0.40, 0.45, and 0.50, the latter being the largest value attainable for all p, see Fig. 6.Yet again, we find good agreement with the analytical results of Blees and Leyte [33] for the case of spheres.As for BCC, the estimate performs better for lower values of p and for lower values of ϕ.Also akin to FCC, the effective diffusivity is monotonically decreasing as a function of p.Note that for ϕ = 0.50, p = ∞ is not included, because the pore phase is discontinuous in that case (a checkerboard structure) resulting in zero effective diffusivity.Neither the simulation nor Eq.
(5) can capture this case.For the SC case, we perform further investigations to find the p value that maximizes effective diffusivity.For the highest solid volume fraction ϕ = π/6, we perform simulations for p = 3, 3.25, 3.5, ..., 9.75, and 10, see Fig. 7.The estimates are noisy, so to estimate the position of the maximum more accurately, we fit the completely  (14) to the data using least squares, providing an excellent fit in the investigated range.This fit indicates that the maximum occurs at p max ≈ 5.4, and is equal to D(p max ) ≈ 0.791.The fact that the two-point bound in Eq. ( 2) is ∼ 0.793 for this solid volume fraction suggests that the class of superballs closely approximates the optimal shape in SC lattices in terms of effective diffusivity.However, it cannot be excluded that some other class of shapes, such as cubes with rounded edges, better approximates the optimal shape.The estimates further suggest that the ζ 2 parameter captures the behaviour around the optimum accurately.That SC lattices have this qualitatively different behaviour from BCC and FCC is not completely surprising considering that SC, as p increases, provides straight, axis-oriented "highways" of pore phase that are very beneficial for transport.For large p approaching a cube, this benefit is partly cancelled out by the superball having too sharp edges.These 'highways' do no appear in the same manner for the other lattice structures.Another relevant property of the lattices is that, as mentioned before, ϕ max is an increasing function of p for SC, but decreasing for BCC and FCC.

V. CONCLUSION AND FURTHER WORK
We have performed simulations of obstructed diffusion in periodic lattice configurations of non-overlapping, impermeable superballs to study effective diffusivity as a function of shape and solid volume fraction for simple cubic (SC), body-centered cubic (BCC), and facecentered cubic (FCC) lattices.In the simulations, point particles diffuse in the intermediate space between the solid superballs.For SC lattices, we find that for constant solid volume fraction, optimal effective diffusivity is obtained for a particular superball with p ≈ 5.4 between a sphere and a cube, whereas for BCC and FCC lattices, no such optimum exists.
Further, the estimate for effective diffusivity using the microstructural parameter ζ 2 captures the qualitative impact of shape in lattice structures accurately, but the quantitative result is sometimes incorrect i.e. for FCC clusters of cubic and near-cubic particles.However, the fact that the trends are captured suggests that the correlation functions contain the relevant information but that the estimate could be refined.One manner in which to do that would be to consider weighted integration of the integrand in Eq. ( 4), and replace ζ 2 in Eq. ( 5) )] β(r, s, θ)dθdsdr, (15) where β(r, s, θ) is a weighting function.Determining the optimal β can be understood as a functional regression problem, solved through least squares fitting and cross-validation.
This requires evaluating the correlation functions on a prescribed grid as opposed to the approach proposed herein, and would be an interesting prospect for further work.Infrastructure for Computing (SNIC).

FIG. 1 .
FIG. 1. (Color online) (SINGLE COLUMN FIGURE, 80 mm) Examples of effective diffusivity curves, showing the D t /D 0 ratio, i.e. the obstruction factor, as a function of time t.Asymptotic values D ∞ /D 0 are obtained as end points of these curves.The curves shown are from simulations in lattices of spheres.
the range [log(r min ), log(r max )].To generate (one set of) sample values of the correlation functions, first three different starting points are sampled uniformly in [0, L] 3 (one for each function).Second, two random unit vectors are sampled, used for determining the directions in which s 2 (r) and s 2 (s) are evaluated.Third, a random set of r, s, and θ is sampled, creating a triangle for evaluation of s 3 .Using a uniformly random unit quaternion converted to a rotation matrix representation, this triangle is rotated randomly.We use lowercase notation to indicate that the correlation functions are treated as binary stochastic variables here.
FIG. 3. (Color online) (SINGLE COLUMN FIGURE, 80 mm) The maximum attainable solid volume fraction ϕ max as a function of p for SC (yellow; bright gray), BCC (orange; medium gray),

FIG. 5 .
FIG. 5. (Color online) (SINGLE COLUMN FIGURE, 80 mm) Effective diffusivity (in units of D ∞ /D 0 ) in BCC lattices for ϕ = 0.20 (upper, yellow; bright gray) and ϕ = 0.25 (lower, orange; medium gray), showing simulated effective diffusivity (filled circles and lines) together with the estimated effective diffusivity (empty circles).Results for spheres (p = 2) are in good agreement with those of Blees and Leyte [33] (empty squares).Note that the shape of the curves are distorted by including p = ∞.

ACKNOWLEDGMENTS
Zheng Ma and Salvatore Torquato are acknowledged for fruitful discussions and comments.The financial support of the Swedish Research Council (grant number 2016-03809) is acknowledged.The computations were in part performed on resources at Chalmers Centre for Computational Science and Engineering (C3SE) provided by the Swedish National

FIG. 7 .
FIG. 7. (Color online) (SINGLE COLUMN FIGURE, 80 mm) Effective diffusivity (in units of D ∞ /D 0 ) in SC lattices for ϕ = π/6, showing simulated effective diffusivity (filled circles) together with the estimated effective diffusivity (empty circles) and the fit of the phenomenological model (solid line).