Numerical estimation of transport properties of cementitious materials using 3 D digital images

A multi-scale characterisation of the transport process within cementitious microstructure possesses a great challenge in terms of modelling and schematization. In this paper a numerical method is proposed to mitigate the resolution problems in numerical methods for calculating effective transport properties of porous materials using 3D digital images. The method up-scales sub-voxel information from the fractional occupancy level of the interface voxels, i.e. voxels containing phaseboundary, to increase the accuracy of the pore schematization and hence the accuracy of the numerical transport calculation as well. The numerical identification of the subvoxels that is associated with their level of occupancy by each phase is obtained by increasing the pre-processing resolution. The proposed method is presented and employed for hydrated cement paste microstructures obtained from Hymostruc, a numerical model for cement hydration and microstructure simulation. The new method significantly reduces computational efforts, is relatively easy to implement, and improves the accuracy of the estimation of the effective transport property.


Introduction
Transport processes in cementitious materials play a crucial role in both the degradation process of building materials as well as in the containment of hazardous wastes.Challenges include difficulties in the overall description of the multiple scales that characterize the microstructure of this reactive porous media, covering a wide range of sizes, i.e. from nanometer-sized pores to centimeter-sized aggregates.Even when only considering a cement paste microstructure, the pore sizes cross at least four orders of magnitude (10 nm -100 µm).When it comes to the description of transport processes in porous media, a parameter that is of paramount importance is the effective transport, e.g.diffusion, coefficient.An efficient and reliable method is needed to obtain the effective diffusion coefficient (D eff ) with an acceptable accuracy.A number of experimental methods already exist but most of them are either time consuming or have limitations.A representative starting point for the numerical analysis of transport properties in cement based microstructures can be obtained by either using: 1) experiments, e.g. by X-ray computed tomography (CT) [1,2], or 2) advanced 3D simulation models.In the numerical hydration models, the resulting microstructure can be simulated as a function of the degree of hydration, particle size distribution, chemistry, water to cement ratio, and reaction temperature.The two main methodological approaches in modeling the microstructural evolution of cement-based systems are: 1) digitalization of the experimental image of the real microstructure, e.g.CEMHYD3D [3], and 2) a particle-based continuous concept of growing spheres, e.g.IPKM [4], Hymostruc [5,6], μic [7].CEMHYD3D simulate hydration and microstructure evolution based on a cellular automata movement and phase change of each discrete voxel.The model is limited to represent the microstructure of hydrated cement paste at levels smaller than the voxel size used, which also affects the ability to describe the pore structure into a small detail.Refining the resolution of these systems is often desirable but this directly increases the simulation computing time significantly.On the other hand, the particle-based models have a vector approach and the size of the developing pore structure is not limited by resolution of the digitalisation.Theoretically, there is no imposing resolution limit on the microstructural formation and pore space representation by growing spheres models [7].However, calculation of effective properties of obtained resolution-free 3D images requires meshing, generally regular lattice voxelisation with a uniform spacing.Numerical modelling, then, is based on a 3D digital image built up from voxels, where each voxel is identified by a phase flag that indicates which phase is represented by a certain voxel volume.The size of the system volume, as well as the resolution, have an exponentially increasing effect on the computation cost for finite element (FEM) or finite difference (FD) numerical transport calculation.Generating a realistic structure for a computational grid, which accurately represents heterogeneities while maintaining the computational efficiency, is still an unsolved problem.
The main aim of the research effort reported in this article was to address a novel multi-scale resolution method that maintains the simplicity of structured grids, while extending the resolution limit and keeping the calculations both more accurate and efficient.The proposed method is tested by means of calculating the effective diffusion coefficient of virtual hydrated cement paste microstructures obtained from Hymostruc [5,6].The method is general in nature that should be applicable to a diverse range of investigated porous medias, e.g.microstructures obtained by different numerical or experimental methods.

Numerical implementation
A finite difference (FD) based program for solving the steady state transport problems is written in the C++ programing language within the Hymostruc platform.Hymostruc´s graphical interface was updated to visualize the transport simulation results.Implementation of the FD transport model without employing pre-processing multi-scale algorithm is described first.A virtual 3D microstructure is discretized into a regular 3D mesh (e.g.Figures 1 and 2 a).Each voxel in a lattice is assigned to be a capillary pore or solid according to the particular position in the microstructure.Then, for each sharing surface between the neighbouring voxels in x, y, and z directions a connectivity coefficient was assigned (c x , c y , and c z , respectively, Figure 2 a) and stored in three c vectors (whose lengths correspond to the number of voxels in the system, N).Six neighbour connection were used, considering that the central node is connected by the sharing faces of a cube in x, y, and z direction.Connectivity coefficients are obtained from microscopic transport properties of the central FD voxel and neighbour FD voxels according to the series connection of two conductors, Eq. 1. 1 / (0.5 0.5 ) where k = 1, w, or (w h) helps to represent a shift of the voxel numbering i to the neighbour voxel in x, y, and z direction, Figure 2 a.The FD second order scheme is used to discretize the Laplace's equation.The Laplace's equation in finite difference form, at each node i is given by Eq. 2.

 
, 1 x i i y i w i w z i wh i wh Assembling all these N equations (Eq.2) for all the FD nodes forms a system of global equations which can be represented in a matrix notation Eq. 3.
 A u b (3) where u is the voltage vector (with size of the total number of voxels in the system, N), A is a sparse and symmetric matrix with 7 diagonals (each voxel has 6 nearest neighborings) that contain information about connectivity coefficients of all the bonds among the voxels, and b is the vector of knowns (i.e.boundary condition at the top of the last layer: u (x,y,z=d+1) = 1, Figure 1 b).The obtained system of equations (Eq. 3) is solved by a conjugate gradient algorithm with an optimized matrixvector multiplication.Next, the flux in z direction at each node i, J i,z was obtained by solving the Fick's first law.The FD scheme used for this is: The effective steady state flux J eff is calculated by averaging the fluxes in z direction in the total system volume (total number of nodes N = w h d).Then D eff can be obtained from the calculated effective flux (normalized to the flux through the same system dimensions without any solid inclusions, J 0 ) according to Eq. ( 5).
The relative diffusivity is the ratio of the effective diffusivity (D eff ) of a diffusing specie in a porous media relative to its value when diffusing in bulk water (D 0 ), and ranges between 0 and 1.

 
, , , 0.5 ( ) ( ) The simulation of transport properties is implemented in two different ways: 1) transport through capillary pores only; or 2) transport through both capillary pores and hydration products (more precisely through the CSH gel).For those capillary porosities that still exceed the capillary pore percolation threshold, effective diffusion is dominated by the capillary pore space because its diffusion coefficient is about 400 higher than the coefficient of CSH gel [3].Non-hydrated cement grains and portlandite (CH) are considered to be impermeable.Transport is influenced by both transports through capillary pores as well as CSH gel pores.This is because of the layered nature of the CSH gel which connects capillary pores.This paper, however, focuses only on the transport through capillary pores.The multi-scale modeling including the transport through CSH gel should be further investigated with the presented multi-scale resolution refinement method.

Multi-scale resolution refinement
In the proposed multi-scale resolution refinement algorithm, first the sub-voxel level information is obtained by pre-procesing refinement of the resolution.This sub-FD-voxel information is then transfered to the FD algorithm (Figure 2).A schematic impression of a fractional quantification of the FD interface voxels is depicted as examples in Figures 2 and 3. First, the pre-procesor for obtaining the sub-voxel information is described.The fractional composition of each of the FD voxels is quantified by increasing the pre-processing resolution (here by up to 83.3nm/voxel).The porosity (i.e.fractional occupancy) within a FD voxel is calculated by counting the pore voxels.The voxel's state (solid or pore) is considered to be represented by the center of the pre-processing voxel.For the example shown in Figure 3 b there are 22 solid voxels so the porosity is P = (64-22)/64.In this approach, the very small particles (as depict in example Figure 2 b), as well as the boundaries of the bigger particles, are now schematized by a much better accuracy.Only the information on the composition of FD voxel is transferred, while the pore morphology (e.g.connectivity in all directions) of the sub-(FD-)voxel level is not accounted for on the sub-voxel level.
Next, passing of the sub-voxel information to FD algorithm is described.In this first approach, the transport property of the interface voxels with the fractional composition is taken to be equal to the overall fractional porosity (e.g. if a voxel has 40% porosity, a relative microscopic diffusivity value of 0.4 is assigned to this voxel).Connectivity coefficients are then calculated from the microscopic transport properties of the neighbour FD voxels according to Eq. 1.Therefore, for interface voxels the connectivity coefficients (c i ) are calculated more accurately by considering their fractional compositions.Moreover, it may be argued that the morphology of interface voxels is also somewhat considered/averaged with a nearest-neighbour voxels (Eq.1); e.g. the connection to the completely solid voxel has finally a zero connectivity (Figures 2 and 3).

Simulation plan
The proposed method is tested for an effective transport calculation of virtual hydrated cement paste microstructures obtained by Hymostruc.Hymostruc simulates the development of the 3D microstructure as a function of the particle size distribution, water-cement ratio, chemistry and temperature.Initial un-hydrated cement particles ranging from 50-1 µm are randomly placed in a cubic volume of 100 μm 3 , which is usually considered as representative for modelling purposes [8].Around 38000 particles were thrown in an envelope, with a particle size distribution following the Rosin-Rammler equation with Blaine value 400 m 2 kg -1 .Simulations are done for Portland cement with a water to cement mass ratio of w/c = 0.4 and 0.3 at 20 °C.During hydration, the outer expansion of the cement particles is calculated according to the so-called particle expansion mechanism [20,23] that is also accounting for the shell overlaps.In order to obtain microstructures with drastically different pore morphologies, two porosities were considered: nominally 11% and 3% (as calculated with the finest resolution of 83.3 nm/pixel).Example of simulated microstructure of the hydrated cement paste is shown in Figure 1 a.Location of solid particles is described by six parameters in 3D Cartesian coordinate system, i.e. the centre coordinates of particles (x, y, z), the diameter of non-hydrated cement grain, and the layer thicknesses of hydration product (Figure 1 a).The simulated microstructure is digitized to form a 3D matrix of cubic voxels with a number of imposed resolutions.Examples of 3D numerical results are presented in Figures 4 a and b.A cubic volume of 100 μm 3 was considered as representative for transport modelling purposes [8].
The employed FD numerical resolutions for the transport modelling were: 0.25, 0.33, 0.5, and 1, that corresponds to a FD system size of 400 3 , 300 3 , 200 3 , and 100 3 voxels, respectively.Furthermore, the effect of pre-processing multi-scale resolution refinement (for obtaining fractional volumes for the boundary voxels and transferring this information to FD property assignment) is also investigated.The finest pre-processing resolution used was 83.3 nm/pixel, resulting in a system with a size of 1200 3 voxels.The digitized 3D structure of the capillary pores was evaluated with the flood fill algorithm to determine its pore connectivity.With this algorithm, for a 3D structure the pore neighbours can be found in a 6 direction configuration or in a 26 directions configuration.In this paper, a six-direction connectivity is used (same as for FD analysis), where the neighbouring voxels share six planes with the central voxel.The degree of capillary pore connectivity is calculated from the number of the connected pore voxels divided by the total number of the pore voxels.
Figures 5, a and b, depict the effect of pre-processing and FD resolution on the calculated relative effective transport property (and relative error) for simulated hydrated cement paste microstructures with two different nominal porosities, P = 11% and P = 3 %, respectively.Each curve corresponds to a different FD resolution used to calculate the effective transport property.The reference value for calculating the relative error was the one obtained with the highest FD resolution (i.e.0.25 μm) and the finest pre-processing resolution (83 nm).The calculated error referenced to this value is shown on the right axis (ordinate) of the Figures 5, a and b.
First, the numerical results obtained for microstructures with the nominal porosity of 11 % are discussed (Figure 5 a).It can be seen that for the FD model predictions calculated without employing the pre-processing algorithm (represented by the right end point of each curve), the errors tend to decline with increasing FD resolution, however, in a somewhat scattered way.With increasing the fineness of the pre-processing resolution up to 1/12 μm (i.e.83 nm), the errors of the predictions tend to decline.Figure 5 a also shows that the calculated relative effective property values do converge.Figure 6 a shows the effect of resolution on the calculated porosity.The calculated porosity is also decreasing with increased resolution.The discrepancy from the decreasing trend of the effective transport with increased resolution can be explained by the change in pore morphology as well as by a possible change in the surface fraction area in the boundary condition with different resolutions.The effect of the resolution on the percentage of connected porosity of these virtual microstructures is shown in Figure 6 b.There is a significant increase in the connectivity of the pore structure with increased resolution.This is as expected because of the more correct description of the fine pores, which can now add to the more connectivity within the network of pore voxels.The pathways that would seem to be depercolated at lower resolution are percolated at higher resolutions when smaller pores can be resolved.
Next, the numerical results obtained for hydrated cement paste microstructures with a nominal porosity of 3 % (Figure 5 b) are discussed.For this low porosity it can be observed that the effective diffusion coefficient obtained without employing the pre-processing algorithm, now tend to increase with increased FD resolution.These results are represented by the right end point on each curve.The increase in effective transport with increased resolution can be related to the remarkable increase in connectivity of the pores in the microstructure.This effect of the resolution on the percentage of connected porosity of virtual microstructures is shown in Figure 6 b.When employing resolutions of 1 and 0.5 μm/pixel, the percolation threshold is reached, and the microstructure depercolates (no open porosity, i.e. 0 % connectivity).However, with increased resolution there is a significant increase in the connectivity of the pores.Figure 6 a shows the effect of resolution on the calculated porosity.The calculated porosity is also decreasing with increased resolution.Therefore, the influence of the pore connectivity on the effective transport at these low porosities now overrules the opposite influence of porosity (compare Figures 6 a and b).By combining these two opposing effects of connectivity and porosity change with increased resolution (Figure 6), while considering the effective transport, one would expect a reverse U shape type of curve.Remarkably, this type of curves are really observed in Figures 5 a and b.With increasing fineness of the pre-processing resolution up to 1/12 μm (i.e.83 nm), the errors of the predictions tend to decline, i.e. the calculated effective property values finally converge.Figure 7 depicts an exponential increase in CPU time with increased resolution.The time for solving the effective transport largely depends on the complexity of the microstructure (compare the results in Figure 7 for aimed microstructures with 11 % and 3 % porosity).The promising finding of Figures 5, a and b, is that it is possible to obtain a high accuracy prediction by using 'coarse' FD resolution, e.g. 2 μm/pixel, if high resolution preprocessing information is incorporated in the transport model.Therefore, a simulation run that would normally take many hours, is now possible in less than a minute, with a similar high accuracy.This demonstrates the powerful capability of the presented multi-scale resolution refinement method.The sub-micro level resolution used in the reported calculations was up to 83.3 nm/pixel, yielding to a system with a size of 1200 3 voxels.The pre-processing in this finest resolution took 150 seconds to run.
The next step in the planned research is to investigate a multi-scale modelling system and also to try to include the morphological influence on the transport properties when passing the information between the modelling scales (nano, micro, meso and macro).  TM) i7 CPU 950 @ 3.07GHz, 6GB RAM).

Conclusion
A multi-scale resolution refinement approach is proposed to mitigate a resolution problems for transport simulations in porous materials.The simulation results showed a dependency of the effective diffusion on a morphology of the investigated structure.Increasing of resolution results in a significant increase of the pore connectivity.Pathways that would seem to be depercolated at lower resolutions turned out to be percolated at higher resolutions, since in this case the smaller pores could also be resolved.On the other hand, the calculated porosity is decreasing with increased resolution.With increasing the fineness of the pre-processing resolution up to 1/12 μm (i.e.83 nm), the values of the predictions first tend to incline, and later start to decline and converge.Combining the opposing effects of connectivity and porosity change with resolution explains the reverse U shape type of curves (Figure 5 a and b) for the effective transport assessments obtained by increasing the pre-processing resolution.Further improvements of the method would be to include the morphological influence of transport properties when transferring information between modelling levels.Further research will focus on an extension of the method towards a multi-scale model that combines nano and micro (cement paste) level, as well as mortar and concrete scales.
a) simulated microstructure (grey-cement, red-inner hydration product, yellow-outer hydration product), b) steady state flux (J) across z-axis, adiabatic boundary conditions employed on 4 side faces parallel to flux.

Fig. 2 .Fig. 3 .
Fig. 2. a) FD implementation, position and size of coordinates: width (x), height (y), and depth (z), b) preprocessing refinement: each sharing surfaces between neighboring FD voxels in x, y, and z directions has an assigned connectivity coefficients c x , c y , and c z , respectively.calculated by Eq. 1.

4 .
3D distribution of molecule concentrations in hydrated cement paste at steady state: a) P nominal = 11 % and b) P nominal = 2.5 %.

5 .
Effect of pre-processing resolution and FD resolution on calculated effective transport property and relative error, for simulated microstructure: a) P nominal =0.11, w/c=0.4,b) P nominal =0.03, w/c=0.3.
CPU time requirement for calculations of effective diffusion property of different virtual microstructures (Intel (R) Core