Modeling and computational homogenization of chloride diffusion in three-phase meso-scale concrete

Abstract A computational homogenization technique for modeling diffusion in concrete is introduced with emphasis on the influence of the aggregate content and variability. The highly heterogeneous material is investigated on different scales by combining Variationally Consistent Homogenization on numerical microstructures with analytical techniques accounting for lower, unresolved, length scales. The concrete structure consists of the cement paste, the embedded aggregates, and the Interfacial Transition Zone (ITZ) in between the two. Diffusion takes place in the cement phase, as well as in the ITZ. Since the thickness of the ITZ is, typically, much smaller than the diameter of the aggregates, the effect of the ITZ can be modelled as a surface transport around the aggregates. The occurrence of different aggregate sizes is described via the Particle Size Distribution for given sieve curves, as described in design codes. The Particle Size Distribution curve is split into two parts. The effect of smaller aggregates is homogenized analytically using a mixture rule. This results in an effective matrix material consisting of cement paste and the smaller aggregates. Synthetic structures are then generated numerically to account for the larger aggregates. At first, a dense sphere packing is created based on the Particle Size Distribution. This information is used to generate a weighted Voronoi diagram, which is modified by a shrinking process. This procedure allows us to create periodic Representative Volume Elements for numerical investigations. The overall diffusivity of the concrete mixture is evaluated upon using Variationally Consistent Homogenization, in the context of Finite Element analysis, for the generated RVEs and compared with analytical homogenization results and experimental data. It is found that, depending on the Particle Size Distribution, the ITZ has a large effect on the effective properties.


Introduction
Concrete is an important material for building structures such as e.g. retaining walls or bridges. It is highly heterogeneous on different scales. Concrete is often used as composite material in form of reinforced concrete to increase e.g. tensile capacity and ductility. Due to its heterogeneity on the small scale and the thereto related porosity, chloride ions can diffuse in the concrete and come in contact with the reinforcement bars causing chloride induced corrosion. Analysis and discussion on the behavior of corroded bars in concrete were presented e.g. by Auyeung et al. [3] and Zandi et al. [28]. Another type of rebar corrosion, namely carbonation induced corrosion, results from the ingress of atmospheric carbon dioxide (CO 2 ) into the concrete. The process of carbonation of concrete and its effect on reinforcement bars are discussed e.g. by Anstice et al. and Hussain et al. [1,11] On the meso-scale, concrete consist of three phases: Cement paste, mineral aggregates and a zone around these having higher porosity and higher diffusion properties. This Interfacial Transition Zone (ITZ) and its influence on the effective diffusivity of concrete was recently examined in experimental and numerical studies. Much of the experimental research has focused on identifying and evaluating the ITZ properties and it's thickness, which is found to lie in the range of 0-50 lm, e.g. [7,27,26]. Analyses of the ITZ under different aspects including the aggregate content and its impact on the effective diffusivity were carried out by Gao et al. [10].
In numerical modeling, computational homogenization is a common tool to account for the multiple length scales of concrete. One specific approach is Variationally consistent Computational Homogenization (VCH), e.g. [12,13], where the homogenization of the transient problem can be carried out in a direct fashion. The basis for numerical investigations is knowledge about the micro-and meso-structure of concrete, i.e. the spatial distribution of the constituents in space. A well established method to account for the aggregate content in concrete is the use of sieve curves. There exists a number of studies which analyse the effect of aggregate content on the overall properties of concrete, see e.g. [4,14,18,21,23]. In those contributions, the meso-scale aggregates are resolved as randomly distributed spherical inclusions. However, the need of resolving the aggregates leads to high computational costs if a large variety of aggregate sizes is considered. Therefore, the before mentioned numerical studies are restricted to rather homogeneous distributions of aggregate sizes, small aggregates are ignored. Furthermore, Caré [5] reported that, in particular the small aggregates are important for the overall diffusivity of chloride ions through meso-scale concrete. Hence, ignoring small aggregates in numerical simulation will lead to results that underestimate the effect of the ITZ on the overall diffusion properties of the compound. Moreover, several studies, e.g. [6,25], reported the importance of the shape of the aggregate which cannot be represented by meso-scale models with spherical aggregates. Therefore, the application of Voronoi diagrams for generation of meso-scale models with randomly distributed aggregates might be a suitable alternative for more natural aggregate shapes as recently studied e.g. for poly-crystalline structures [9] and for meso-scale asphalt concrete [22,24]. Finally, it was reported in several studies that the size of the ITZ depends on the aggregate size [16,17,19].
In this paper, we aim at executing efficient numerical multiscale modeling of chloride diffusion in concrete to gain insight about the corrosion process in reinforced concrete. We generate artificial structures for numerical simulation using Particle Size Distributions (PSD) for sieve curves obtained from the literature [5,26]. To overcome the aforementioned deficiencies of previous studies we take into account a large variety of aggregate sizes as well as irregular aggregate shapes. The generation of structures undergoes three steps and is discussed in more detail in Section 2: First, the generation of a dense sphere packing is performed by means of the Lubachevsky-Stillinger algorithm [8]. Second, a weighted Voronoi tessellation [20] is generated. Third, the (initially convex) polyhedral cells stemming from the Voronoi tessellation are shrunk in a random fashion to mimic the irregular shape of mineral aggregates, e.g. crushed rock. The synthetically generated aggregate structures are used to simulate mass diffusion within a control volume. The transport of ions in the ITZ is modeled as diffusion along the surface of the mineral aggregates whilst the aggregates themselves are assumed to be impermeable for ion diffusion. The upscaling of the diffusion properties from the meso-scale control volume towards the overall macroscopic diffusivity is executed in terms of VCH in Section 4. To account for the ITZ effect associated with small aggregates, the aggregate content is split into smaller and larger aggregates. Whilst the larger aggregates are geometrically resolved via Voronoi tessellation, the diffusion on the scale of the smaller aggregates is homogenized analytically via appropriate mixture rules in Section 3, i.e. without resolving the smaller aggregates geometrically for the sake of reduced computation times. In Section 5, the multi-scale model is tested and validated against experimental data taken from the literature.

Preliminaries
We model mass diffusion in three-phase concrete by applying VCH. In this approach, the effective properties, such as the diffusivity of chloride through the three-phase concrete, are obtained from homogenization over a Representative Volume Element (RVE) 1 . In this Section, we propose a technique to generate such RVEs artificially for numerical simulation. The generation procedure is based on sieve curves taken from the literature. Such RVEs have to be, on the one hand, small enough to allow for numerical simulations at reasonable costs. On the other hand, they have to be large enough to be representative for the full structure. For the investigation of three-phase concrete, we assume separation of scales between the macro-level and the RVE level, subsequently called meso-level. Thus, L ) L Ã ) L M , see Fig. 1.
Moreover, we introduce an even smaller length scale, denoted micro-scale, that consists of the smaller aggregates, ITZ and cement paste. By contrast, the matrix material on the meso-level represents a mixture stemming from homogenization of the micro-scale properties. Our motivation for introducing this a priori split into larger aggregates (within the mesoscopic RVE) and smaller aggregates (within a "micro-RVE") is to ensure reasonably low computational costs. More precisely: For computational efficiency, the number of aggregates that can be fully resolved in the 3D RVE is limited. If the range of aggregate sizes that needs to be resolved in one RVE becomes too large, a huge amount of smaller aggregates is needed to match a desired sieve curve. Decomposing the sets into smaller and larger aggregates helps to bypass this limitation and to achieve suitably simplified yet realistic aggregate structures.
Our goal is to implement a numerical procedure that is able to generate synthetic aggregate structures based on given sieve curve data. An example sieve curve is shown in Fig. 2 a) with aggregate sizes ranging between d = 0:075 mm and d = 4:75 mm [26] 2 Here, the black dots and lines represent upper and lower bounds for the grading of a concrete mixture from ASTM Type I Portland cement, sand and water. Their volume/mass are summed up and presented according to the cumulative mass fraction given in m-% versus the mesh size that is directly related to the average aggregate size. The sieve curve of the synthetic structure shall range within those bounds. As mentioned above, the explicit modeling of all aggregates would lead to unacceptably complex RVEs which are not easily accessible for numerical simulations at reasonable costs. We, therefore, homogenize the material properties on the micro-scale a priori and, for convenience, divide the sieve curve at an amount passing of 50 m-%. Only the aggregates above the corresponding mesh size are resolved in the mesoscopic RVE. Hereby, the diameter of the synthetic aggregates is computed as the diameter of the volumeequivalent sphere. The sieve curve of an example realization of syn-1 Note that, throughout this paper, quantities referring to properties of the mesoscopic RVE will be labeled with the subscript "Ã". 2 Note that this data set, which is chosen for illustration purposes, corresponds to mortar rather than to meso-scale concrete. However, it is important to remark that the method introduced in this paper is general and also suitable for much larger aggregates.
thetic aggregates is given in Fig. 2 a). The (geometrically resolved) large aggregates are marked with red dots. The (unresolved) small aggregates for an example micro-RVE problem are given as blue dots. Note that, for the synthetic structure, each dot corresponds to one single aggregate whilst, in a standard sieve curve, it is counted how many aggregates can pass through that particular mesh. The aggregates are impermeable and, in consequence, not modeled explicitly but taken as cutouts of the RVE as shown in Fig. 2 b). One advantage of the artificial generation of structures is that the RVEs can be a priori generated in a periodic fashion to facilitate the application of Periodic Boundary Conditions (PBC) during upscaling. The ITZ is assumed to be very thin and its thickness to be much smaller than the aggregate length scale. Thus, we technically describe diffusion through the ITZ on the aggregates surface and do not resolve the ITZ domain explicitly.
The volume fractions of the concrete phases are defined as where LA; SA and M denote the large or small aggregates and the homogenized matrix material containing the pure cement paste and the small aggregates. The volume of the RVE is jX Ã j, and the subscript A refers to all aggregates (SA and LA). The volume fraction of the ITZ is assumed to be much smaller than n A ; n M . It is important to remark that, even though the volume fraction of the ITZ is negligibly small, the ITZ has a strong effect on the overall diffusivity of the structure. Due to the choice of diffusivity of the ITZ we do not specify the thickness but assume that this thickness is very small and in the range of lm.  The mesoscopic RVEs are generated in three steps that are explained in more detail subsequently. At first, a dense sphere packing with given PSD is generated. Second, a weighted Voronoi tessellation is generated to model the artificial aggregates. Third, these Voronoi cells are shrunk to adjust the desired volume fraction.

Generation of dense sphere packings and polyhedral aggregates
We start by generating dense periodic sphere packings as precursor structures for the desired mesoscopic part of the sieve curve. We use the Lubachevsky-Stillinger algorithm [8,15] to find a dense packing of spheres with variable radii. The algorithm calculates moving elastic spheres with random initial positions and velocities that grow under a given PSD. The algorithm uses event-driven molecular dynamics shown in Fig. 3 a). At time t ¼ 0, each sphere i is located at position x i , has an initial velocity v i and an initial diameter d i ¼ 0. Based on the given sieve curve, we specify a PSD and take it as input for the growth rate g i of the spheres.
The spheres grow linearly in time with Þand may undergo elastic collisions. Here, t n and t nþ1 are the times at which collision n and n þ 1 occur. If the spheres i and j get in contact, their velocities are recomputed. The time t nþ1 is calculated from Hereby, no solution means no collision, one solution means a contact between the spheres, and finding two solutions means that there is an overlap of the spheres resulting in a re-computation of velocities. The process is finished if the time step between two collisions is smaller than a predefined value. An example for the resulting periodic dense packing is presented in Fig. 3b).
In a second step, a radical Voronoi tessellation [2,20] is performed on the basis of the sphere packing. The positions of the spheres and the radii serve as input parameters for the generation of the Voronoi diagram. Thus, the generator points, i.e. the center points of the spheres, are assigned a weight according to their radius [22]. The resulting geometry has an aggregate volume fraction of n A ¼ 1 and therefore needs a change to adjust the desired volume fraction.

Adjustment of the volume fractions
The last step is about shrinking the Voronoi cells. The Voronoi diagram contains the three dimensional Voronoi cells in form of polyhedrons and the associated seed points P S;i . To adjust the volume fraction, we re-use the spheres stemming from the dense sphere packing in the previous step and randomly distribute N k points on this sphere. Thus, we create points on a spherical hull around the seed points. The Voronoi cell and the hull may intersect each other, see Fig. 4a). The distribution function of these random hull-points P H;k is uniform. We use spherical coordinates in the form where r k ¼ 0:5 d k is given by the PSD, h k is the latitude, and u k is the azimuth angle. If a point of this hull P H;k , k = 1; 2; . . . ; N k , is located outside the related Voronoi cell, the point is moved towards the center of the sphere. This process is controlled by the shrinking factor s 2 0; 1 ½ . Starting from s ¼ 1, this factor is successively reduced in an iterative fashion until all hull points are located within the associated Voronoi cell and until the volume fraction of aggregates in the RVE matches the desired value. In other words, if a hull-point P H;k is located outside the associated Voronoi cell, the radius of this point is multiplied with the shrinking factor via r i;k ¼ s r i . Here, the index i refers to the ith Voronoi cell. Subsequently the related coordinates x k are used to recalculate the volume of the new shrunk (and possibly non-convex) polyhedron. That means, each hull point is multiplied by the same shrinking factor, but not all points are shrunk. The randomness of this procedure allows to generate aggregates with highly variable shapes as illustrated in Fig. 4  Note that this procedure is a priori performed in a periodic manner within a unit cube, see Fig. 4 b). Since we consider diffusion through the volume of the cement paste as well as via the ITZ, i.e. the aggregate surface, the size of the RVE in relation to the aggregate diameters becomes relevant. It is, therefore, necessary to scale the artificial microstructure. This is effectuated based on the largest sieve size d max and the largest generated sphere with diameter d unit max . The volume of the RVE is computed as 3. Meso-scale model of chloride diffusion

Preliminaries
We execute homogenization of the micro-scale based on a Taylor assumption to capture the complexity of the structure of concrete. As mentioned in Section 2, the aggregates below a specific mesh size d min = 0:2 mm are homogenized analytically on the micro-scale to derive material properties for the numerical study on the meso-scale. The microscopic control volume X M contains the cement paste as matrix material with a distribution m R ð Þ of spherical small aggregates with radius R i per volume mixture. Each sphere is surrounded by the ITZ. Using the Taylor assumption, we seek the overall diffusivity D M in cm 2 /s from averaging the flux due to a uniform gradient f = $c. We therefore need to define at first the diffusion coefficients for all phases.

Bulk and interface diffusion
Chloride diffusion is strongly dependent on the diffusivity of the three phases: aggregates, cement paste, and ITZ. We assume the aggregates to be impermeable, i.e. D SA ¼ 0 cm 2 /s. The diffusivity of the cement paste is set to D CP = 1EÀ7 cm 2 /s which is in the range of realistic diffusion properties according to [5,18]. The diffusivity of the ITZ is higher than the one of the cement paste. We introduce the scale factor f to easily assess the diffusivity of the ITZ and define the ''effective diffusivity"D I;S aŝ Note that the factor f (dimension of a length) is a material parameter that contains, in an implicit fashion, information about the thickness of the ITZ as well as about the intrinsic diffusivity of the ITZ. Whilst the thickness of the ITZ could be detected e.g. via 3D imaging, the (inhomogeneous) distribution of the intrinsic diffusivity within the ITZ is very hard to assess. Therefore, we propose to use f as the only model parameter for the effective ITZ diffusivity. Subsequently, we will demonstrate how f and, thereby, the effective diffusivity of the ITZ can be identified from experimental data. Furthermore note that Nilenius et al. [18] proposed the comparable formulation b D I;S = D I;S t = D CP f , where D I;S is the intrinsic diffusivity and t characterizes the corresponding thickness of the ITZ whereby f has dimension of a length. The latter formulation bears the advantage that f is the only unknown material parameter in the model instead of the generally two unknown parameters t and D I;S .
In order to determine the effective diffusivity, we need to derive the volume fractions of the constituents. The volume fraction of the cement paste contained in X M;Ã isñ CP ¼ 1 Àñ SA , in which the volume of the ITZ is negligible and the volume fraction of small aggregates is defined as Here, the volume of one small aggregate V SA;i ¼ 4=3 pR 3 i is computed as the volume of a sphere with radius R i . 3 This volume fraction of the small aggregates in X M;Ã is linked to the volume fraction of small aggregates in X Ã , see (1), viañ SA ¼ n SA =n M .
We introducen Ã I;S in cm À1 as a weight factor that quantifies how much the ITZ, surrounding a spherical aggregate, contributes to the overall diffusivity. Note thatn Ã I;S is introduced to take into account that, given a concentration gradient f c ½ , only those parts of the sphere surface contribute that are parallel to f c ½ . Those parts are called "active parts" of the surface. The derivation ofn Ã I;S is given subsequently.

Homogenization of the micro-scale properties
The relations for quasi-static mass diffusion conditions on the micro-scale are given as whereÎ = I À n n is the surface projection with the normal vector n. The material parameters are given in Table 1, and the diffusivity tensors are written in the case of isotropy as We adopt the assumption of uniform gradient, i.e. f c ½ ¼ f M in X M;Ã to derive the effective diffusivity using the Taylor assumption where the distribution m R ð Þ is the frequency of aggregates per volume in the mixture jX M j; V R ð Þ is the volume of a sphere (representing a small aggregate) andñ CP ¼ 1 Àñ SA . We derive the volume of one small aggregate as the volume of a sphere V R ð Þ ¼ 4 Using the surface projection, we identify the ''active" partŜ R ð Þ of the surface aŝ Hence, we only evaluate that part of the aggregate surface, i.e. of the ITZ, that is parallel to the uniform concentration gradient f M . Considering spherical aggregates, we restrict to the simplest Thus, we can derive the weight factor as The frequency m R ð Þ depends on the PSD, see Fig. 2a). It is calculated based on the sieve curve, i.e. each spherical aggregate is unique with radius R i , which gives for any function F R ð Þ, where N is the number of all aggregates on the micro-scale. This leads tõ Assuming isotropy, we rewrite the effective mass flux upon combining (6) For the subsequent numerical studies, we introduce the normalized effective diffusivity as which serves as input parameter for the meso-scale RVE problem. In (21), we set D SA ¼ 0.

Strong and weak format of the meso-scale problem
We develop a model for mesoscopic mass diffusion in a threephase concrete under static conditions. We assume the large aggregates to be impermeable, i.e. D LA ¼ 0, and compute the mass 3 Remember that the small aggregates are not resolved geometrically. Instead, a Taylor rule is applied which assumes the small aggregates to be spheres.  ð22cÞ c I ¼ c p I on L I;D ; ð22eÞ where n is defined as the outwards pointing normal of the aggregate surface. In (22d), J c M ½ Án is the leak-off from C I into X M . We assume continuous concentration which leads to the interface condition To simplify notation, we henceforth skip the indices and use c to denote the concentration in X M Â X A and on C I . The constitutive relations for the mass flux J in g/(cm 2 s) are defined as where square brackets r ½ relate to operational dependency of on r, and f T c ½ is the gradient of the mass concentration c. We apply the gradient operator f c ½ for the matrix material and the tangential gradient operator f T c ½ :¼ I À n n ½ Áf c ½ that is used on C I . The diffusivity tensors D in cm 2 /s andD in cm 3 /s are, in case of isotropy, defined as where, again, b D I = f D CP . We introduce the standard space-variational format for the system in (22) as follows: Find c x ð Þ in the appropriately defined spaces C M Â C I that solve where C 0 M and C 0 I are the appropriately defined test spaces.

First-order homogenization in the spatial domain
We are now able to define the two-scale problem based on the single-scale problem in Section 3.4. This problem depends on running averages in the weak format and scale separation via firstorder homogenization which are defined subsequently. At first, we replace the space-variational problem in (26) by that of finding c x ð Þ 2 C FE 2 that solves where the pertinent space-variational forms in (27) are given as These forms represent running averages on cubic domains X Ã with the side length L Ã which are located at each macro-scale spatial point x 2 X. The underlying volume averaging operators are defined as h i Ã;M : where the volume averaging operator over X LA vanishes as the aggregates are impermeable, i.e. D LA = 0. The weak form in (27) is the basis for the VCH discussed subsequently. We introduce scale separation and first-order homogenization. Thus, we decompose the sub-scale mass concentration c into a macro-scale part c M and a fluctuation part c s within each RVE in the format Note that, in this stationary problem, the macroscopic concentration c is energy-free and plays the role of a ''rigid body" mode. Therefore, c ¼ 0 can be chosen for convenience. The two-scale trial and test space are defined as where C and C s Ã;i are the finite element discretized trial spaces for the macro-and meso-scopic problem on the RVE derived as i.e. there exist c s with corresponding space C s Ã;i for each single RVE We restate the two-scale problem taking into account the interpretation dc s i ¼ dc s , by finding c; c s

Macro-scale problem
We test now with a purely macroscopic test functions, i.e.
As a result, we obtain the homogenized macro-scale problem from (33) as where the macro-scale flux is defined as

Meso-scale problem on an RVE
Taking into account the difference operator s t Ã x ð Þ ¼ x þ ð ÞÀ x À ð Þ, we introduce the additional model assumption that c s is periodic which makes the RVE problem uniquely solvable. Here, x 2 C þ Ã are the image points and x 2 C À Ã the corresponding mirror points on opposite parts of the RVE surface. We apply periodicity in a variational format such that where the RVE forms are defined as Here, k M and k I represent boundary fluxes on C þ Ã;M and L þ I . The test space for boundary fluxes is defined as with the space L 2 of square integrable functions and T Ã;I represents k i 2 R for each print L þ Ã;I . By this, we state the space-variational RVE format as follows: Solving (40) allows us to derive the macroscopic diffusivity tensor from the (linear) relation such that the components of the macroscopic diffusivity tensor in a Cartesian system can be, due to the linearity of the problem, computed as Note that, in general, the macroscopic diffusivity can be anisotropic. For convenience, we define the normalized and isotropic macroscopic diffusivity as where D ii is the trace of the macroscopic diffusivity.

Analytical homogenization of the meso-scale problem
For comparison reasons, we subsequently define alternative formulations for D norm that are based on analytical mixture rules of the meso-scale properties. First, we adopt (21) for the meso-scale problem such that for D LA ¼ 0, where the surface weight factor is defined aŝ Here, N is the number of (spherical) aggregates in the RVE. Note that the definition (45) is computed in a way that only those parts of the aggregate surface, i.e. of the ITZ, that are parallel to the overall concentration gradient contribute to the diffusion process.
A second alternative is to compute the macroscopic diffusivity from the mixture rule Here, the surface weight factor is defined aŝ Again, N is the number of (spherical) aggregates in the RVE. The surface weightn I is computed in a way that the entire surface of the aggregates, not only the "active" part, is contributing to the diffusion process.

Numerical investigations
For our numerical investigations, we generate various random meso-structures according to the sieve curves given in [5,26]. All material parameters are specified in Table 1. They are based on experimental data taken from literature [5,18,26]. We execute VCH of the meso-scale problem with geometrically resolved large aggregates whilst the smaller aggregates are only considered via the mixture rule (21). The reason for this split into (resolved) large and (unresolved) small particles is that in typical sieve curves, e.g. that one given in Fig. 2, the largest aggregates are more than one order of magnitude larger in diameter than the smallest aggregates. Fully resolving the entire bandwidth of small and large aggregates would make it extremely costly and demanding to resolve the full structure and to solve the RVE problem. The numerical results of this hybrid technique (micro-problem -analytical homogenization, meso-problem -VCH) are validated against the experimental data. Moreover, we compare the VCH solutions with results obtained from applying analytical homogenization rules for both, micro-and meso-scale problem according to (44) and (46). Subsequently, we investigate the properties of our multi-scale model for diffusion in three-phase concrete. In Section 5.1, we study how sensitive the predicted macroscopic diffusivity is for variations in the sieve curve. In Section 5.2, we validate our model against experimental data from the literature. Note that, for convenience, we compute the diameter d of a mineral aggregate as the diameter of the sphere with the equivalent volume.

Model sensitivity
We investigate the sensitivity of the multi-scale model for changes in the underlying sieve curves. We consider four different scenarios based on the sieve curve data included in Fig. 5. The bounds for the sieve curves correspond to those in Fig. 2 a). In the sieve curves, each blue dot represents a small aggregate and each red dot a large aggregate in a random example realization. Hereby, discrete sieve curves have been used for both, analytical homogenization and VCH. The resulting effective diffusivity is plotted versus the volume fraction of mineral aggregates for three different realizations per data point which are equipped with an error bar to account for the scatter of the results. We vary the scaling factor from f ! 0 to f ¼ 0:05 mm.
First of all, we observe in all result curves the general trend that, for f ! 0, the macroscopic diffusivity decays more or less linearly with the volume fraction of the aggregates. This behavior is due to the fact that the mineral aggregates are assumed to be nonconductive (D SA ¼ D LA ¼ 0). Hence, the overall diffusivity becomes smaller the larger the aggregates' volume fraction is.
In Fig. 5a), the structure under investigation only contains small aggregates, i.e. the meso-problem is homogeneous, and only the mixture rule (21) is used. Aggregate diameters are chosen in the range d 2 0:075; 0:2 ½ mm. We observe a strong increase of the overall diffusivity if f is increased even though the overall value of f ¼ 0:05 mm is very small compared to those that will be used in Section 5.2. For the studied random realizations of the microstructure, the error bars are negligibly small.
Vice versa, the structure analyzed in Fig. 5 b) only contains large aggregates with d 2 0:3; 1:18 ½ mm diameter. The overall diffusivity is studied in terms of VCH and the mesoscopic mixture rules (44) and (46). We observe that, for all chosen values f, the mixture rule (46), followed by the mixture rule (44), predicts the highest overall  diffusivity. Moreover, we observe that the error bars for the results obtained from the mixture rules are negligibly small. This is a consequence of the fact that no meso-structural information, neither aggregate shape nor spatial arrangement of aggregates in the RVE, are taken into account in this case. By contrast, spatial resolution of the large aggregates leads to a significant scatter in the results as seen in the error bars for the VCH results at f ! 0. This scatter indicates that the RVE size, chosen for the sake of numerical efficiency in this analysis, is not truly representative, see the discussion in Section 2.1. In this case, the aggregate shape and the spatial distribution of the aggregates in the RVE influence the computed overall diffusivity significantly, in particular at higher aggregate volume fractions. By contrast, an increase in f, i.e. a setting where diffusion through the ITZ is the dominating mechanism, reduces the sensitivity of the overall diffusivity for the distribution and shape of individual aggregates in the RVE. Thus, the macroscopic diffusivity for the studied random realizations is almost identical. But it is most important to remark that, overall, the increase in diffusivity for increasing f is much less pronounced than for the micro-structure investigated in Fig. 5 a).
This can be observed in Fig. 5 c) in similar fashion, where the previous micro-structures (d 2 0:075; 0:2 ½ mm) and mesostructures (d 2 0:3; 1:18 ½ mm) are combined. Here, the overall diffusivity is augmented in comparison to b) but less than for the purely microscopic structure in a). In other words, the behavior in c) is a combination of a) and b). Finally, Fig. 5 d) shows results for an approximately mono-disperse meso-structure (d % 1 mm). The absence of smaller particles on the micro-and on the mesoscale leads to a situation where the diffusion via the aggregate surfaces, i.e. via the ITZ, is of minor importance for the studied values of f.
To illustrate the effect of the ITZ on the meso-scale diffusion pattern we plot the normalized mass flux, i.e. the diffusive ''speed" of the ions, in Fig. 6 a) in the absence of any ITZ (f ! 0) and in b) for f ¼ 0:5 mm. We observe that, if the ITZ is absent, the diffusion mechanism tries to circumvent the particles. This leads to strong local concentration gradients in the cement paste, mirrored in a heterogeneous mass flux with strongly curved stream lines. By contrast, the ITZ is the preferable zone for diffusion. Therefore, the mass flux profile in Fig. 6 b) is rather homogeneous, and the stream lines are attracted by the ITZ.
Altogether, this first study shows that the overall diffusivity of three-phase concrete is driven by mainly two concurring effects: First, high volume fractions of non-conducting mineral aggregates lead to lower diffusivity. Second, the smaller the particles, i.e. the larger the surface-to-volume ratio, the larger the overall diffusivity of the three-phase concrete. Clearly, the latter is a length scale effect. We conclude that the ITZ surrounding the small aggregates is carrying the most significant part of the diffusivity.

Model validation
In the second study, we aim at validating our multi-scale model against experimental data from the literature. The question to solve is: Can we find, for a given D CP , a suitable value for f such that the numerically predicted overall diffusivity matches the experimental data? For that purpose, we extract experimental results from [26,27,5] as shown in Fig. 7.
Hereby, the experimental results by Yang et al. [27] are obtained for a concrete mixture with the sieve curve displayed in Fig. 2. Taking into account the spreading of the experimental data, a very good approximation of the experimental results can be obtained by choosing a very small value f ! 0. In other words, the surface diffusion effect associated with the ITZ in the mixture used in [27] is of minor importance. However, the approximately linear drop of overall diffusivity with increasing aggregate volume fraction in properly matched by the multi-scale model.
In [5], three different mortar mixes with varying aggregate contents are studied. The three mixtures are denoted ''fine" (d 2 0:315; 1 ½ mm), ''medium" (d 2 0:315; 2 ½ mm) and ''coarse" (d 2 2; 4 ½ mm). Further details of the sieve curve are not reported. Therefore, we assume a sieve curve similar to that of [27], but scaled to similar diameter ranges as the ones called fine, medium and coarse. This scaling is effectuated according to (5) and leads to the ranges indicated in Fig. 8.
The procedure to compare the numerical multi-scale model to the experimental results is as follows: We use the available data for the coarse mixture to adjust the parameter f of our model. We find that the choice f ¼ 0:5 mm provides a good agreement of the VCH model with the experimental data for the coarse mixture. The same parameter f ¼ 0:5 mm is now used to predict the overall diffusivity for the medium and the fine mixture. We find that for the medium mixture, the VCH prediction is matching the experimental data very well. For the fine structure, VCH results in a slight underestimation of the diffusivity of approx 10%.
Applying the same parameter f ¼ 0:5 mm for the mixture rule (44), we find a mild overestimation of the overall diffusivity. By contrast, the mixture rule (46) leads to a significant overestimation of the data and is, for the sake of lucidity, not displayed.
In view of the incomplete data on the sieve curves used in the literature, the agreement of both methods, VCH and analytical homogenization only taking into account the active parts of the  ITZ, i.e. those being parallel to the macroscopic concentration gradient, is remarkably good. It requires further experimental data to decide whether it is necessary to carry out the numerical effort for VCH or whether it is sufficient to evaluate the (much cheaper) analytical homogenization rule.

Conclusions and outlook
The aim of the current study is to establish a modeling technique for chloride ion diffusion in three-phase concrete consisting of smaller and larger mineral aggregates embedded in a cement paste matrix with the ITZ. The ITZ is a thin layer of increased porosity and diffusivity surrounding the aggregates. We propose to use a hybrid multi-scale approach based on sequential analytical homogenization of a micro-scale problem (containing small aggregates and cement paste) and computational homogenization of the meso-scale problem (containing large aggregates and the homogenized matrix material from the micro-scale problem). To be able to deal with realistic microstructures, we present a method that allows to use sieve curve data and to generate, in an automated fashion, random (periodic) RVEs representing the meso-scale of three-phase concrete. The algorithm is based on (i) a generation of a dense sphere packing according to the required PSD, (ii) a weighted Voronoi tessellation, and (iii) a shrinking procedure of the Voronoi cells to further randomize the structure and to adjust the volume fraction of the mineral aggregates.
Having these tools at hand, we are able to analyze the overall diffusivity of three-phase concrete for a large variety of sieve curve data. We, therefore, first investigate how sensitive the multi-scale model is for changes in the sieve curve. The main finding is that the overall material behavior of three-phase concrete shows a clear size effect in the presence of the ITZ: The smaller the aggregate size the larger becomes the overall diffusivity. At the same time, large volume fractions of (impermeable) aggregates can reduce the overall diffusivity significantly. The question which of these concurring effects is dominant is decided by the thickness and the diffusivity of the ITZ.
In a second study, we validate the numerical multi-scale model against experimental data from literature. Hereby, calibrating the model for one given data set allows us to predict the behavior of mixtures with different aggregate content with reasonable accuracy. However, more complete experimental data is required for a more detailed validation.
Altogether, the proposed multi-scale modeling approach gives an insight into the diffusive chloride ingress in three-phase concrete. The relevance of the ITZ is clearly supported by the current findings. As to ongoing and future work, we plan to successively address the subsequent aspects: carry out time-dependent (transient) investigations, investigate the effect of irregular aggregate shapes, extend the structure generation algorithm towards dependence of ITZ thickness on aggregate size Furthermore, we aim at investigating if the proposed method can be used to inversely identify the intrinsic diffusivity of the ITZ if geometrical data about ITZ thickness is available, e.g. from X-ray Computed Tomography. Finally, it could be interesting to study a type of percolation effect that occurs when the ITZs surrounding two close aggregates start to interact.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.