Discrete superstructures in low-resolution images Open Access

Computer simulations can serve as powerful tools to understand and predict real-life results only when they accurately mimic experiments. In particular, simulations of fluid flow in complex geometries are highly relevant to a variety of scientific and industrial applications.^1–14 Computer simulations originating from first principles (based on the solutions of partial differential equations (PDEs) such as the Navier–Stokes equation for fluid flow) can be used to verify laboratory results and guide experiments.¹⁵ The complexity of real-life problems requires the solution of PDEs to be numerical, which in most cases necessitates the use of discrete meshes. The basic uniform Cartesian mesh is not only involved in the solutions of PDEs, but is also routinely found in the imaging of real-life objects in the form of pixels or voxels. For example, two-dimensional digital photography produces a set of pixels while three-dimensional (3D) X-ray, magnetic, or optical scanning produces a set of voxels.

The finite memory of digital computers limits the number of mesh elements (e.g., pixels) in an image, and it seems that an optimal approach to image an arbitrary-shaped object is to distribute the mesh elements in space uniformly along Cartesian directions. This approach generates pixels in two dimensions or voxels in three dimensions. However, mapping real-life objects or processes onto a Cartesian mesh unavoidably leads to the discretization error. In the context of flow simulations, the discretization error impacts both the geometry representation and the flow field, with the latter originating from the numerical solution of PDEs. As the result, this discretization error contaminates the physics of flow simulations. Minimization of the discretization error is required to validate the computer models, and ensure they accurately reproduce the experiment. Refining the mesh (increasing the resolution) is an option; however, this route leads to prohibitive scanning and/or simulation times. Therefore, obtaining a smaller discretization error at lower resolutions is of utmost practical importance.

Previous research on pore-scale flow simulations has reported^7,16–18 improving flow simulation accuracy via (i) the addition of detailed information about each pore (i.e., a solid boundary between mesh nodes), and/or (ii) pores of the sample to be sufficiently resolved by the mesh. Low-resolution images with unresolved pores are naturally avoided for accurate flow simulations. Low-resolution computed tomography (CT)-scanned images are also impacted by the fundamental limitation of image contrast, which further complicates the entire process of obtaining accurate geometry followed by flow simulations.^15,19 In this study, we deal with the packings of closely packed spheres discretized on a uniform cubic mesh. We refer to the discretization resolutions of below $∼10$ voxels per sphere diameter as “low,” while resolutions of above $∼50$ voxels per sphere diameter are considered as “high.” If we consider the hydraulic diameter $dh$ to be about one third of a sphere diameter for a packing with porosity of 0.35,^20,21 then the low resolutions will correspond to $dh≲3$ voxels while high resolutions to $dh≳15$⁠.

Here, we analyze and aim to minimize the discretization error in simulations of flow through a porous sample, providing the new physical insights into flow simulations. The pore-scale flow simulations are highly sensitive to the topology (flow “pathways”) and morphology (pore “shapes”) of the pore space. In this study, we (1) apply fractional discretization and visualize superstructures, which can be visible when regular sphere packings are mapped onto the discrete mesh at low resolutions, (2) simulate flow to obtain permeability of a porous sample, (3) establish similarity in the permeability error between regular and irregular geometries, (4) vary the free adjustable parameter of the lattice Boltzmann flow solver, which controls the discretization errors (“magic number”), (5) highlight the existence of the “null point,” where low-resolution flow field provides accurate permeability value due to self-cancelation of the discretization error contributions, and (6) propose a linear correlation between the geometrical parameter $χ$ and the magic number. Finally, we demonstrate the presented linear correlation leads to an excellent error reduction of the permeabilities computed from the three-dimensional (3D) images of laboratory-prepared samples relative to the experimentally measured values.¹⁵ 

Conventional pore-scale flow simulations require discretization of a porous sample by mapping it onto a cubic mesh and marking each mesh voxel as either solid or fluid based on its center location relative to the sample solid phase [Fig. 1(b)]. The ratio (number of fluid voxels)/(total number of voxels) defines the discrete porosity. We initially define the discretization resolution as the number of voxels per sphere diameter. For regular geometries and low discretization resolutions, we maintain the discrete porosity of each geometry close to its analytical value with minor adjustments of the sphere radii during discretization, if there is a noticeable difference between analytical and discrete porosities.

To visualize the discrete meshes originating from regular geometries, we consider the smallest case of a porous sample to be a single sphere. When coupled with periodic boundary conditions, this geometry results in a simple cubic (SC) packing with an analytical porosity of $1−π/6≈0.476$ [Fig. 1(a), top row]. The selected flow simulation approach requires an integer number of nodes (L) per each dimension of the simulation domain. We use the periodicity property and replicate the SC unit cell U times along each Cartesian direction. When the ratio $L/U$ is non-integer, each cubic unit cell has a non-integer number of mesh nodes per edge, while L always remains an integer. We refer to this process as the fractional discretization procedure.²² The top row of Fig. 1(b) illustrates that the fractional discretization procedure results in feature-rich voxelized geometries (⁠ $U=2$⁠, $U=3$⁠) compared to the geometries of similar resolution with an integer $L/U$ ratio (⁠ $U=1$⁠). Note that increasing U alone adds no new information to the analytical geometry—it is a simple replication of SC unit cell. An identical outcome is observed for both body- and face-centered cubic packings (BCC and FCC), as seen in the middle and bottom rows in Fig. 1. These results show that the fractional discretization results in superstructures—structures with dimensions that significantly exceed one unit cell.

To assess the accuracy of flow simulations, we focus on the permeability, calculated using the average velocity in the direction of the applied pressure gradient. (Permeability quantifies the capacity of a given geometry with voids to conduct a fluid.) The simulated flow is single-phase and pressure-driven, and occurs in the voids of geometries formed by the closely packed spheres at various void space fractions (porosities). We simulate a zero-Reynolds number flow which obeys Stokes PDE, with solutions obtained using the two-relaxation-time lattice Boltzmann method (LBM)²³ implemented as described in Ref. 24. The no-slip boundary condition is enforced using the bounce-back rule. Applied pressure drop and the corresponding macroscopic flow are directed along one of the principal axes. We perform flow simulations in the void space of SC and irregular (containing 14 400 spheres, Fig. 2 inset) geometries of identical porosities, and calculate the permeability error relative to a reference value. Both SC and irregular geometries have their own $∼0.1$%-accurate reference permeability values $kref$ obtained using extrapolation,²⁴ see Tables I and II for the complete list of values. These reference values can be obtained in different ways, for example pushing resolutions to prohibitive levels such as $∼103$– $104$ voxels/sphere diameter, replacing bounce back with higher-order boundary conditions and using moderate resolutions of $∼102$⁠, or using any non-LBM numerical scheme which will provide resolution-free permeability values based on the solution of Stokes equation.

Figure 2 compares variation of the relative error in permeability vs resolution for the SC and irregular geometries. This figure demonstrates that increasing the number of unit cells per domain edge U, while preserving a non-integer $L/U$ ratio, reduces the scatter of the relative error in permeability for the SC geometry. Also, with the increase of U, the relative error for the SC geometry begins to follow the irregular one. This means that when a superstructure within the SC geometry reaches a sufficient size, the SC geometry displays the resolution–permeability error dependency similar to the irregular geometry. This finding suggests the existence of superstructures in not only regular but also in irregular geometries.

To visualize the impact of superstructures on the flow field for the SC geometry, we color each fluid voxel according to its absolute velocity magnitude. To reveal the skeleton of the flow field, Fig. 3 shows about 200 voxels with the highest magnitude. The skeleton in Fig. 3 resembles the features of the superstructures shown in Fig. 1(b). Figure 3 reveals that (1) similar to the superstructures seen in Fig. 1(b), the size of each flow field skeleton also significantly exceeds one unit cell, and (2) the skeleton (and superstructures) do not resemble the pore space of the underlying analytical geometry. The first point suggests that construction of any scheme to numerically solve Stokes PDE inside a unit cell cannot capture the entire superstructure because of insufficient information on the pore space. For the considered flow problem, we are not aware of any numerical scheme constructed outside of a unit cell: the numerical schemes such as the lattice Boltzmann and finite difference rely on the information from a given voxel plus its direct neighbors, which is below the unit cell scale. The second point demonstrates that the superstructures—not pores—control flow at low resolutions, and implies that varying the resolution changes both morphology and topology of the voxelized pore space. Conversely, with increasing resolution from low to high values, the impact of superstructures disappear and flow occurs through the pores of the underlying analytical geometry (see Fig. 8 as an example for SC geometry).

We note that visualization of the flow skeleton and the corresponding superstructures for irregular geometries will be limited as any observed local velocity maxima (which form the superstructures in Fig. 3) can be attributed to a slightly larger pore sampled by a given voxel. However, similarities in the resolution–error curves in Fig. 2 suggest that irregular geometries also contain the superstructures.

The discretization error is the key artifact separating computer simulations from their real-life counterparts. A free parameter known as the “magic number,” $Λ$⁠, controls the spatial discretization error in two-relaxation-times lattice Boltzmann simulations.^23,25 In this section, we provide essential background details on $Λ$⁠. In later sections of this study, $Λ$ together with geometrical parameter $χ$ will be used to construct a universal correlation, which significantly reduces the discretization error in permeability simulations.

LBM simulates the fluid with fictitious particles that occupy the discrete mesh and propagate along the prescribed discrete links at discrete time steps. On each iteration, the particles collide at mesh nodes according to a predefined collision operator. The collision operator can be formulated differently,^26–28 but it always includes at least one collision (relaxation) rate. In the basic case of Stokes flow and the collision operator with a single rate, this rate controls both viscosity of the simulated fluid and the spatial discretization error.²⁶ The adjustment of viscosity separately from the discretization error can be done with at least two collision rates, which resulted in the formulation of two-relaxation-times (TRT) collision operator.²³ Then, $Λ$ is a specific combination of the LBM collision rates.²⁵ Note that not all collision operators with two rates separate the viscosity control from the discretization error (see Sec. 2.1.3 in Ref. 22). For the collision operators with multiple rates,²⁷ several combinations of the collision rates need to be fixed to separate viscosity adjustment from the spatial discretization errors; see discussion in Secs. 2.1.4, 2.2.3 in Ref. 22.

The collision operator defines the numerical scheme in the bulk, away from the solid–fluid interface, while the boundary condition completes the scheme at the interface. The continuum description of flow around a solid obstacle assumes zero-flow velocity at the solid boundary which is also known as the no-slip boundary condition. To implement this boundary condition at the solid–fluid interface in voxelized images, we use the popular “bounce-back” LBM boundary condition.^29, $Λ$ can be seen as the parameter controlling (i) the location of the zero-flow boundary (and the corresponding pore “width”) between solid and fluid voxels,²² or (ii) the magnitude of the discretization error or its contributions (see Fig. 3 in Ref. 24; also Eq. (15) in Ref. 30, where parameter $τ$ can be interpreted as the one impacting the discretization error).

The existence of the free parameter $Λ$ requires it to be assigned a value before running a simulation, as there is no clear guidance for a particular choice of $Λ$ for simulations in complex geometries. The numerical permeability obtained with LBM and the bounce-back rule for a given geometry is arbitrary, and it is controlled by $Λ$⁠,^22,25 although the impact of $Λ$ on the permeability decreases with the mesh refinement. If we consider popular reformulations of the collision operator, such as the Bhatnagar–Gross–Krook operator-based²⁶ (BGK), multiple-relaxation-times,^27,31 or cumulant-based²⁸ [e.g., Eq. (12) in Ref. 32], the free choice of $Λ$ impacts all of them. That is, the choice of $Λ$ is of fundamental importance to obtain accurate simulation results.

Currently, a robust theoretical analysis for the simple system of flow between two parallel plates suggests taking $Λ=3/16$ (or $3/8$⁠) for the exact flow field at any discretization resolution in a horizontal (or $45°$-inclined) channel relative to the underlying mesh.^22,33 Similarly, $Λ=1/8$ provides the exact average velocity, canceling the velocity integration error.^34, Figure 4(a) adds the impact of $Λ$ to the results from Fig. 2, showing the error in permeability vs resolution at different values of $Λ$⁠. Figure 4(a) also shows that once the geometry includes curved boundaries, these $Λ$ values no longer result in the most accurate permeability.²² Note that in Fig. 4(a) increase in the resolution from low to high result in all curves converging to zero from above, crossing the zero error value, and then slowly continuing to converge up from below. [This also includes $Λ=3/8$ in Fig. 4(a), as can be seen in Figs. 8(c) and 8(d) in Ref. 22] This is counterintuitive because increases in resolution may result in larger permeability errors.

The impact of $Λ$ on simulated permeability can be significantly reduced^18,22 after replacing the first-order bounce-back boundary condition with a higher-order one. This replacement is only possible after adding information on the location of the pore boundaries between the mesh nodes to the LBM simulations. However, this approach has fundamental limitations when recovering pore surfaces from real voxelized images. Namely, (i) recovery of a surface using, e.g., the marching cubes algorithm³⁵ does not converge to the analytical result even in the case of a single sphere;³⁶ and (ii) during imaging, contrast loss occurs non-uniformly within each pore,¹⁵ which prevents pore surface recovery from the intensity of each voxel.

Each curve in Fig. 4(a) demonstrates a distinct point which we call the “null point” or P0. At P0, the resolution can be low (unresolved geometry and flow field), the error in permeability is small (⁠ $∼1$%), and any resolution deviation from P0 increases the error magnitude. This small error at P0 originates from the self-cancelation of three components of the discretization error: the LBM scheme away from the solid–fluid interface, the boundary condition, and the integration while calculating the average flow rate. Figure 4(a) also shows that P0 is different for each $Λ$ value, while it is similar for the SC and irregular geometries. P0 is present but not discussed in other studies.^{7,22,34,37,38}

Please note that this linear correlation is obtained for the three basic regular packings: FCC, BCC, and SC. For completeness, we also add all periodic packings with the porosities of 0.2 and 0.15. Here, at lower resolutions we see a deviation of FCC packing from the general trend. These two porosities are omitted for the linear fit in Fig. 4(b).

We assess the accuracy of the linear correlation given by Eq. (4) for irregular geometries using six irregular computer-generated packings and two laboratory-packed CT-scanned samples. This accuracy check is based on the parameter $χ$ obtained for each voxelized geometry and running LBM flow simulations with $Λ*$ calculated according to Eq. (4).

Computer-generated geometries include periodic irregular packings of mono- and poly-dispersed spheres as well as a packing confined laterally by the wall of a cylindrical container. The confining wall imposes partial ordering of sphere locations and introduces porosity and flow velocity maldistribution near the wall, propagating 3–5 diameters from the wall into the bulk (see Figs. 2 and 4 in Ref. 40) The impact of the confining wall is significant for the ratio of 10 sphere diameters per cylindrical container, and therefore, this geometry is also used to assess Eq. (4).

We also employed experimental geometries from our previous study¹⁵ for evaluation of Eq. (4). We packed two different types of commercially available glass beads, 0.47 mm beads1 and 0.54 mm beads2, in 8.98 and 9.04 mm glass tubes in water under ultrasonic vibration (Appendix A in Ref. 15) resulting in samples P3 and P4, and determined their permeability experimentally. Hereafter, both samples were scanned using CoreTOM X-Ray CT scanner (XRE Tescan, Ghent, Belgium) using the tube voltage of 60 kV and power of 15 W. Each sample was scanned at 18 resolutions. The number of two-dimensional (2D) projections ranged from 250 to 1800, while the exposure to obtain each projection was about 4.5 s, Appendix B1 in Ref. 15. Gray CT images were segmented using global thresholds each equal to the laboratory-determined porosity values. Simulated permeability for samples P3 and P4 agrees with experimental values within 1%. The reference dimensionless permeability (porosity) value for P3 is $5.82×10−4$ (0.354) while for P4 is $6.02×10−4$ (0.355), see Table II and Appendixes A3 and C3 in Ref. 15. More details on the preparation and experimental measurements for samples P3 and P4, their imaging, image processing, and simulations are provided in Ref. 15. All 3D images,⁴¹ 2D CT projections,⁴² and 2D optical scans of beads1 and beads2 with experimental logs⁴³ are available online.

Sphere size distributions (SSD) obtained from 2D optical scans of beads1 and beads2 were used to computer-generate irregular periodic packings at porosity of 0.362 (Fig. 24 in Ref. 15). Extrapolating from the LBM simulations with $Λ=0.05$⁠,²⁴ we determined the reference 0.1%-accurate permeability values for computer-generated packings as $6.851×10−4$ for beads1 SSD and $6.833×10−4$ for beads2 SSD, see Table II.

To benchmark the accuracy of Eq. (4) on CT images of real objects, we need to consider the impact of image contrast. X-ray computed-tomography scanning produces images with a reduced contrast [blue “3” panel in Fig. 5(c)] in comparison with computer-generated and discretized geometries, which have no contrast loss. Recent research¹⁵ has demonstrated that at lower resolutions the contrast of CT images affects the simulated permeability. This was observed using an operator-independent global segmentation procedure of gray CT images, based on the laboratory-measured porosity. The images with maximum contrast [e.g., green “3” in Fig. 5(c)] were obtained by downsampling high-resolution CT images [similar to “64” in Fig. 5(c)], see Appendix B2 in Ref. 15. The same study observed that for these maximum-contrast images the error in permeability is identical to the computer-generated geometries, Figs. 2, 27(a), and 27(b) in Ref. 15. Therefore, we use the images of maximum contrast from the earlier study¹⁵ to simulate flow with $Λ*$ from Eq. (4). These downsampled, maximum contrast images are available online.⁴¹ 

Figures 5(a) and 5(b) show a comparison of the simulated permeabilities obtained in the current study using a common value of $Λ=0.25$ (equivalent to $τ=1$ for BGK collision operator) and $Λ*$ from Eq. (4). For all geometries, the proposed correlation (4) brings the permeability error to unexpectedly low levels, as seen in Figs. 5(a) and 5(b). Our results demonstrate that selection of $Λ*$ according to Eq. (4) enables accurate simulation results from unresolved, highly-voxelized images.

Decrease in the discretization resolution in flow simulations allows to save computational time significantly: the computational complexity of the employed LBM simulation approach scales as $O(resolution5)$ meaning that computational time between resolutions 4 and 64 differ $O(165)≈O(106)$ times. Despite more efficient finite difference solvers for Stokes PDE are available⁴⁴ (although at a cost of a higher discretization errors due to the reduced voxel connectivity), the computational efforts will still grow rapidly with the resolution increase. Also, reducing the resolution saves CT scanning time: according to our experience, scanning a sample at the resolution of 4 voxels per sphere diameter takes about 10 min while the resolution of 64 voxels requires 1200 min. Alternatively, reduction of the CT resolution allows to scan larger volumes of the sample at a fixed scanning time.

To date, efforts to improve the simulation accuracy of flow through porous media have targeted localization of pore surfaces of original, non-discretized geometries.^18,32,45 During imaging of a real porous medium, scanning equipment maps the porous medium geometry onto a discrete uniform mesh. Superposition of the real geometry and the mesh results in the formation of a superstructure (Figs. 1 and 3). The existence of this superstructure suggests that relying on pore surfaces (within one unit cell or a representative volume) is fundamentally limited with decreases in resolution to low values. This limitation originates from the incomplete information about the pore space geometry and the corresponding flow field at the unit cell level. At each low resolution, the superstructure corresponding to this resolution and geometry controls flow. More generally, superstructures control not only flow but also morphological and topological information about a given geometry. Superstructures also prevent pore-level analysis from producing the optimal magic number, $Λ$⁠. By contrast, increases in resolution from low to high values result in both the recovery of pores and the flow through the image of the original, non-discretized geometry (Fig. 8).

Currently, performing accurate simulations of pore-scale flows demands an accurate, visually appealing representation of the pore space [red “64” panel in Fig. 5(c)]. Reduction of the mesh resolution to a few voxels per pore results in a pore space that is highly voxelated and unattractive to the human eye [blue, red “3” panels in Fig. 5(c)]. However, we show that flow simulations on these highly-voxelated images can accurately reproduce the experimental permeability (i.e., the flow physics) because of the existence of superstructures, which retain information about the pore space over scales significantly exceeding a single pore dimension or representative volume of a porous sample.

The geometrical origin of the superstructures indicates that the presented findings are not limited to a particular numerical method (here LBM, see also Fig. 7 for the finite difference flow fields) or to a particular flow problem (Stokes flow). As mentioned previously,²⁵ variation of $Λ$ in LBM simulations is similar to, e.g., a finite difference method, in which the derivative coefficients can be adjusted. Such adjustment changes the order of convergence of the method, as well as its error magnitude. Here one should keep in mind that for the final simulation results the error magnitude is of key importance rather than its convergence rate.

The correlation (4) identifies the free magic parameter $Λ$ such that the individual error contributions (bulk and boundary errors of the numerical scheme and integration error in calculating the average flow rate) of potentially opposite signs cancel each other out. This formal approach can be applied to other phenomena simulated by numerical solutions of PDEs. As we see it, the key ingredients to successfully implement our approach are: a discrete uniform mesh, numerical scheme (here LBM), existence of the null point where the error is zero, and a target integral quantity (here permeability). We expect that the error cancelation can be achieved not only by the spatial replication of a target object (here a unit cell) to obtain a superstructure, but also by the object's behavior in time (e.g., when the object or its features move). Thus, the presented approach can be used for better positioning of arrays of sensors, improved temporal measurements, or better quantification of pixelized images provided, for example, by drones or satellites.

We are grateful for the allocation of computational resources by the Supercomputing Laboratory at King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia. We are thankful to Gabrielle E. Abelskamp for her editorial support. We thank the anonymous reviewers of JFM-2024-1620 submission for their helpful comments. This work was supported through funding to the Ali Al-Naimi Petroleum Engineering Research Center (ANPERC) at KAUST. Contributions: S. K. suggested the existence of discrete superstructures, introduced the concept of null point, wrote manuscript draft, and prepared figures. T. P. organized and supported this research, supervised it and edited the paper. S. K. introduced the fractional discretization in the joint study with Dr. Irina Ginzburg (10.1016/j.jcp.2014.10.038) after I. G.'s suggestion to preserve symmetry of the analytical geometry during discretization.

The authors have no conflicts to disclose.

Siarhei Khirevich: Conceptualization (lead); Supervision (equal). Tadeusz W. Patzek: Funding acquisition (lead); Supervision (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.