Towards stochastic generation of 3D pore network models of building materials

Pore-scale-based prediction of the hygric properties of porous building materials is on the rise as an attractive alternative for the current experimental procedure. Pore-scale simulations do however require a complete pore network model for the building material. With the currently available characterization techniques, such complete pore network model cannot be established, instead typically fragmented direct (pores sizes, shapes, positions, connections, ...) or indirect (pore size distribution, pore surface area, ...) information is obtained. The aim of this paper is to present stochastic pore network generation, wherein the fragmented pore structure information is used to generate a complete pore network for the building material involved. The novelty of our approach lies in the generation of a PNM by matching the distributions of direct parameters as well as indirect parameters of the input data and the model. Additionally, the position of the pores are no longer bound to a cubic lattice. This workflow will first be tested on a single scale material with a relatively straightforward pore space such as sintered glass. Finally, the hygric properties of the generated network will be compared to the measured properties of the real material as a validation step.


Introduction
Flow and transport phenomena in porous media play a significant role in various fields of science and technology, comprising a spectrum from medical sciences over material sciences to soil and rock sciences. In order to determine the moisture behavior of building components, numerical simulation models are commonly used. However, these models require an good description of the moisture retention and moisture permeability functions, as these are crucial input parameters. The experimental characterization of these hygric properties is impaired by serious weaknesses [1]. The procedure merges ad-and desorption, static and dynamic methods, thus ignoring hysteresis, air entrapment and dynamic effects. The measured outcomes additionally remain incomplete, as they provide not much information about the mid-saturation ranges. Moreover, it necessitates weeks of experiments, and of course requires to have the actual material at hand. This multifactorial group of flaws implies that hygric analysis of existing materials is cumbersome and time consuming.
In order to overcome these deficiencies, a new approach is currently being introduced in building physics. Pore network models (PNMs) have already proven to be useful to determine fluid flow properties in other research fields, such as geological applications [2], fuel cells [3] and drying technology [4], to determine fluid flow properties. PNMs try to represent the studied pore structure as accurately as possible while retaining a certain simplicity by representing the pore space as a network of pore bodies and pore throats. Hence, the PNM tries to capture local features of the pore-space which are important for the fluid storage and transport processes under investigation. Subsequently these PNMs are subjected to invasion algorithms that replicate different (de)saturation procedures: absorption, desorption, imbibition and drying. For unsaturated moisture storage and transport in building materials, Islahuddin and Janssen, developed a hygric pore-scale simulator comprising the coexisting liquid and vapour phases of water [5,6].
Unfortunately, several materials, including building materials, often have a broad spectrum of pore sizes, ranging from nanometers to millimeters. The presence of multiple pore-scales in the studied sample can severely influence its storage and transport properties. Therefore, the applied PNMs need to incorporate information over different length scales. PNMs have the advantage that they are inherently scale invariant, i.e. this technique can be applied to any length interval for which the pore structure has been experimentally observed and analyzed. However, a general approach to combine different PNM at different spatial scales into one global network is currently missing.
As a first step towards multi-scale PNMs, the primary objective of this paper targets the definition of a workflow to generate PNMs based on information gathered from direct 3D imaging techniques such as micro-CT on a single scale. As an example the suggested methodology is tested on a material with a relatively straightforward pore space, a glass filter made by sintering glass particles. Based on 3D micro-CT information the stochastical information about the different PNM parameters such as pore and throat volume, radius and shape descriptor as well as their connectivity is collected. This information is subsequently used to generate a new representation of the pore network, taking into account the correlations between these parameters.
In order to validate the results of the simulation, crucial parameters of the simulation are compared with the measured ones. Additionally the hygric properties of the simulated material are calculated and compared to test the accuracy of the synthetic pore networks.

Methodology
The stochastically generated networks recreate the pore space of the studied material based on input data such as distributions of pore size, throat size, connectivity. These distributions can be obtained from an analysis of pore-space images or in the future even from indirect data acquisition techniques: e.g. mercury intrusion porosimetry or nuclear magnetic resonance testing. The generated networks can be arbitrarily large and hence are not limited by the size of the original image or network.
In the past several workflows have been suggested to generated such network models. They typically progressively build on each other, as they all have the same fundamental structure but each research group tried to incorporate an additional level of complexity to improve the overall simulations. Ioannidis et al. [7] extracted pore and throat size distributions and pore connectivity information from binary images of thin-sections. Starting with regular cubic lattices, they removed nodes and bonds in order to match the desired average coordination numbers corresponding to the simulated porous media. Cubic pores and rectangular throats were then distributed on the remaining nodes and bonds. Subsequently, Arns et al. [8] offered a better method of generating networks, by matching the coordination number distribution instead of just the average coordination number. Finally, Idowu et al. [9] developed tools to preserve the morphological properties of the network. They hence no longer require a cubic lattice, by assessing the following three observations: the minimum and maximum distances between two pores are decided by the user, each throat radius should be smaller than the pores it connects and two adjacent pores cannot overlap. They also allowed to incorporate a first degree of correlation between the building blocks of the network by letting the pore size determine the throat sizes that are connected to them.
In our approach all fundamental information about the pore network such as the inscribed radius, volume and shape factor of pore bodies and throats as well as the coordination number of the pore bodies, is gathered and a statistical distribution is fitted on each of these parameters. Copulas, functions that describe dependencies among variables, are used to this end, as they provide a way to create distributions for correlated multivariate data. Using a copula, a multivariate distribution can be defined by specifying marginal univariate distributions and choosing a particular copula to define a correlation structure between the distributions. Bivariate distributions, as well as distributions in higher dimensions, are possible. Because of their flexibility these functions have become popular for data analysis in recent years.
The entire simulation algorithm can be described by the following steps: 1) Assign the desired dimensions to the random network 2) Place an equivalent number of pores as the original material within this volume.
This can be done randomly or by using a Sobol set of coordinates. 3) Check if the minimum and maximum distances defined by the user are respected 4) Assign a correlated set of geometrical information and coordination number to each pore. Check if no overlap of pore volumes occurs.

5) Define connections between pores by starting with the largest one and connect it
with the i nearest pores, with i equal to the coordination number. For each of the connecting pores, the remaining number of connections is decreased by one. 6) To each connection a weight is assigned based on the radius of the two pore bodies it connects. The smallest throat radius along with other correlated geometrical information is assigned to the branch with the smallest weight.

Results
To illustrate the suggested pore network generation algorithm a sample of ROBU P100 is scanned at a resolution of 2.5 µm. The simplicity of the pore network is confirmed by a narrow pore volume distribution with a diameter between 20 and 100 µm. Moreover the measured open porosity, by using a vacuum saturation setup, of 35.8 % corresponds well with the bulk density of the glass, indicating the absence of closed porosity in the glass particles. Hence, a micro-CT scan with a resolution of 2.5 µm captures the entire range of pore sizes. Subsequently the images are segmented and the pore space is transformed into a PNM using the maximal ball approach developed by Blunt et al. [10]. In this dataset the following independent parameters are present for pore bodies: volume, inscribed radius and shape factor as well as the number of connections with neighbouring pores (coordination number). For the throats, the volume, radius and shape factor are calculated as well as each throat length. The histograms of pore volumes and the throat radii as well as their fitted lognormal distributions are shown in Fig 1A&B. In Table 1 the correlation values for each of these parameters are gathered. Additionally, in Fig 1C, the correlation between the pore radius and pore volume, with a correlation factor of 0.82, is illustrated.  In this example two methods to determine the location of the pore bodies of the PNM are used. In the first method the locations are chosen randomly with the only constraint that two pore bodies cannot overlap. For the second method, the positions are determined based on Sobol sequences, which are quasi random low discrepancy sequences. These sequences use a base of two to form successively finer partitions of the defined interval and subsequently reorder the coordinates in each dimension. This results in a more space filling covering of the simulated volume. In Fig 1D the relationship between the pore radius and connecting throat radii are plotted. Because of the imposed condition that throat radii should always be smaller than the radii of the two pores it connects, all points are situated bellow the diagonal. A good match between the point clouds representing the input data and the output data is achieved.
In Figure 1E the cumulative distribution of the throat lengths are shown for 50 repetitions of both approaches. The red curves, representing the use of Sobol sequences, display less variance, and are also closer to the blue input data curve, compared to the results obtained by randomly choosing the pore positions.
In Fig 1F histograms of the resulting overall porosity value of the simulations based on random and Sobol positioning are illustrated. It is clear that an underestimation of the porosity occurs in the random case. This can be explained by the random positioning, which places the pores somewhat closer together, resulting in more small pores, which have even smaller throat radii and hence a smaller total pore volume. Finally, three simulated volumes with different porosity values are subjected to invasion algorithms as developed by Islahuddin and Janssen, (2017). The moisture storage and transport property curves are calculated for the entire water saturation range (Fig 1 G&H). For both parameters a good match between the original and the simulated pore network is obtained. In Fig 1G the  volumes is related to their respective open porosity. In general all the curves of the simulated volumes respect the curve of the original input volume. Additionally, the saturated permeability value of the simulation is compared with lab measurements. The simulated and measured saturated permeability values are in the same order of magnitude, but the simulated one is slightly higher: 2.12 10 -5 vs. 1.84 10 -5 kg/m s Pa respectively.

Conclusion
In this paper a workflow to generate PNM based on correlated statistical distributions is presented. Two approaches to determine the position of the simulated pores are tested and compared. The results indicate the advantages of using Sobol sequences compared to random positioning. The resulting moisture retention and transport property simulations show a good match between simulations and real measurements for this uncomplicated single scale material.
Because of the broad applicability of PNMs and their versatility to incorporate different scales without exploding the calculation time required for simulating different physical properties, PNMs have the potential to resolve the problems related to simulate fluid flow in complex pore-space systems, such as building materials. In order to incorporate information of different data acquisition techniques such as direct and indirect imaging techniques, as well as the synthesis of different pore networks obtained at different spatial scales (e.g. the integration of nm, µm and mm scale pore networks). This objective will form the bread and butter of our further research, as it allows to combine all available information about the pore space and come up with a full scale network which is necessary to obtain accurate simulations.