A Numerical Study of ITZ Percolation in Polyphase Concrete Systems Considering the Synergetic Effect of Aggregate Shape- and Size-Diversities

The percolation of the interfacial transition zone (ITZ) is generally regarded as an important factor that may accelerate the penetration of aggressive agents in concrete materials, and its threshold is largely determined by the features of aggregates. In most numerical studies about ITZ percolation, both fine aggregates and coarse aggregates are assumed to be the particles of uniform shape, and their size distributions are generally strung together by a single function, which is quite different from reality. To quantify the ITZ percolation associated with the polydispersity of aggregate shapes and size gradations in a more realistic way, the two-dimensional (2D) meso-scale model of concrete is generated by simplifying coarse aggregates and fine aggregates as polygons and ovals, respectively. Moreover, the size gradations of them are also represented by two separate expressions. By combining these models with percolation theory, the percolation of ITZ in the 2D case is explicitly simulated, and the influence of aggregate shape- and size-diversities on the critical threshold ϕagg,c is studied in detail. Based on the simulated results of ϕagg,c, an empirically analytical expression is further proposed to fast predict the ITZ percolation, and its reliability is verified. The results show that the ITZ thickness, average aggregate fineness, coarse aggregate shape, and fine aggregate shapes are the four main contributing factors to the ITZ percolation. Compared with the existing literature, the proposed model here has a broader range of applications (e.g., mortar, concrete, and other granular systems) in the 2D case and can provide the larger predicted results, which may be closer to reality.


Introduction
The interfacial transition zone (ITZ) [1][2][3] is commonly deemed to be an important component around the aggregate particles in concrete, and its distinct characters (e.g., higher porosity, lower strength, and larger permeability) can facilitate the penetration of aggressive ions into materials and the stress concentration to some extent, which in turn leads to the degradation of material performance [4,5]. Specially, the ITZ percolation [6,7] is mainly used to reflect the global connectivity of the ITZs around the aggregates, and its critical threshold can be generally regarded as a contributing factor that determines the micro-structure and macro-properties of materials. According to the percolation theory [8,9], once the global connectivity of ITZs emerges, a dramatic and rapid change in material properties may occur. Therefore, exploring the ITZ percolation in concrete is of great significance for evaluating the effective performance of materials. micro-structure and macro-properties of materials. According to the percolation theory [8,9], once the global connectivity of ITZs emerges, a dramatic and rapid change in material properties may occur. Therefore, exploring the ITZ percolation in concrete is of great significance for evaluating the effective performance of materials.
Due to the complexity and changefulness of the ITZ phase in concrete, it is generally difficult to experimentally identify and quantify the partial characteristics of the ITZ phase, such as ITZ thickness, relative content, percolation behavior, etc. To remedy this deficiency, computer simulation has become a powerful alternative tool in the field of percolation [10][11][12]. Up to now, a large number of researches about the ITZ percolation in concrete have been conducted by simulation, and considerable achievements have been successively obtained. For instance, by assuming the aggregate particles as regularshaped particles (e.g., spheres [13][14][15], ellipsoids [16,17], polygons [18], polyhedrons [19], etc.) with polydisperse sizes, the sensitivity of the ITZ percolation to the aggregate shapes and sizes was studied individually, and the findings were concluded as below: (a) the ITZ thickness is the dominant factor that determines the ITZ percolation [14][15][16][17][18][19]; (b) the percolation threshold of ITZ is very sensitive to the specific surface areas (SSA) of particles and shows an obviously downward trend with the increase of SSA [17][18][19]; (c) with the increase of the aspect ratio of elongated ellipsoids [16,17] or the decrease of the number of polygon edges [18] and polyhedron faces [19], the percolation threshold of ITZ around them could be decreased significantly. All of these studies [14][15][16][17][18][19] supported the importance of aggregate shape and size in the prediction of ITZ percolation. However, there are still several problems that need to be improved because of two simplistic assumptions generally used in the previous papers. The first assumption is that both coarse aggregates and fine aggregates are commonly simplified as regular particles of uniform shape, which overlooks the geometric differences between them [15,16,18,19]. As shown in Figure 1, the aggregate particles in concrete can be roughly classified into two categories: coarse aggregates (i.e., gravels) and fine aggregates (i.e., sands). Morphologically, the geometric shapes of coarse aggregates are very similar to the regular/irregular convex polygons or polyhedrons [18][19][20][21][22][23]. By contrast, the fine aggregates look more like the non-spherical particles of smooth surfaces such as spheres [13][14][15], ellipsoids [16,17], and ovoids [7,24,25], etc. In addition, the other simplistic assumption is that both the size distributions of coarse aggregates and fine aggregates in the literature about ITZ features are generally strung together and quantitatively represented by a single function [15,16,[18][19][20][21][22][23][24], which is also quite different from the reality in concrete. The oversimplification of aggregate shapes and size distributions in previous studies may result in a greater deviation between the predicted approximation and the reality.  To emphasize the effect of the shape diversity of granular components in materials, a series of complex granular systems have been applied recently in the study of ITZ percolation. For example, a combination body of ellipsoids and spheres is adopted to represent the fine aggregates and air-voids in the air-entrained mortar, and the dependence of ITZ percolation on the shape features of both aggregates and air-voids is obtained [26]. All the specific surface areas of the fine aggregate and air void, ITZ thickness, and airvoid volume fraction can be regarded as crucial parameters to the ITZ percolation. In the literature [27], the fine aggregates and lightweight (LWA) aggregates are respectively assumed to be the spherical particles with the fixed thickness of ITZ and the spherical particles with no ITZ, and the sensitivity of ITZ percolation to the relative content of different types of aggregates is presented. The figures in this paper show that the volume fraction of percolated ITZ for the LWA mortar is generally lower than that for the normal mortar. In addition, by assuming the aggregates and fibers as ellipsoids and cylinders, respectively, the effect of the features of aggregates and fibers (e.g., aggregate fraction, gradation, fiber length, etc.) in high performance concrete (HPC) on the percolation of ITZ around them is investigated [28], and the results show that although the ITZ regions in HPC are thinner and not percolated, they can be re-percolated by the addition of just a few fibers. According to the above studies [26][27][28] for the granular system composed of poly-shaped particles here, it can be found that the variation of ITZ percolation in them would be much more complicated than in the mono-shaped particle systems, which also indirectly proves the necessity of distinguishing the geometric difference between coarse aggregates and fine aggregates in the numerical prediction of ITZ percolation in the concrete systems.
In order to overcome these two simplistic assumptions mentioned above and determine the ITZ percolation associated with the polydispersity of aggregate shapes and size gradations, the coarse aggregates and fine aggregates in concrete are respectively assumed to be polygons and ovals based on their real morphologies in Figure 1. Moreover, the size gradations of them are also expressed by two separate functions instead of a single contiguous function as described in refs. [15,16,[18][19][20][21][22][23][24]. On the basis of above two preconditions, the sensitivity of ITZ percolation to the features of both coarse aggregates and fine aggregates in a more realistic way is studied in this paper. Besides, it should be stressed that to reduce the research difficulty and shorten the computation time, the whole work here is conducted in the two-dimensional (2D) case.
The details of this paper are arranged as follows: first, the mathematical expressions of both ovals and regular polygons, the modeling of uniform ITZ around them, and the generating procedure of 2D meso-scale models of concrete are successively presented in Section 2. By combining these models with continuum percolation theory, the effective percolation of the uniform ITZ around both oval fine aggregates and polygonal coarse aggregates is explicitly simulated in Section 3, and the effect of the polydispersity of aggregate shapes and sizes (e.g., the aggregate shape, gradation, and sand ratio) on the critical threshold φ agg,c is quantitatively evaluated in detail then. Based on the simulated results of φ agg,c , an empirically analytical expression is further proposed by the linear regression analysis, and the comparisons of the predicted value here and the numerical results of the existing function are then presented. Finally, some interesting conclusions are highlighted in Section 4.

Geometric Expressions of Ovals and Polygons
The oval is a catch-all term for a family of 2D non-centrosymmetric particles that include circles, ellipses, egg-shaped particles, and other convex-shaped particles with smooth contours, as shown in Figure 2. It has been widely used to characterize the essential features of components in materials [29,30] due to its extensive representation. Based on the geometric features of the oval particles in Figure 2, the profile of ovals can be simply treated as an extension of ellipses, and their mathematical expression can be described by Equation (1). where ζ is a tapering coefficient that varies from 0.0 to 1.0, a and b are the semi-axis lengths of the oval in the xand y-axis directions, respectively, and the ratio of b to a is defined as the aspect ratio of the oval. on the geometric features of the oval particles in Figure 2, the profile of ovals can be simply treated as an extension of ellipses, and their mathematical expression can be described by Equation (1).
The geometric shapes of ovals with different aspect ratio b/a and tapering coefficient ζ.
where ζ is a tapering coefficient that varies from 0.0 to 1.0, a and b are the semi-axis lengths of the oval in the x-and y-axis directions, respectively, and the ratio of b to a is defined as the aspect ratio of the oval. As shown in Figure 3, a set of regular polygons with different numbers of edges are graphically displayed. According to the study in [18], the geometric topographies of these polygons can be mathematically expressed via the collection of the subsets of vertices and edges. Take the regular pentagon as an example; it is assumed that the centroid of the particle is located at the origin of rectangle coordinates and the x-axis passes through one of its vertices. The size and length of the pentagon are denoted by the notation l. Furthermore, all the vertices Pi (i.e., P1, P2, P3, P4, and P5) of a particle can be symbolized by Equation (2). As shown in Figure 3, a set of regular polygons with different numbers of edges are graphically displayed. According to the study in [18], the geometric topographies of these polygons can be mathematically expressed via the collection of the subsets of vertices and edges. Take the regular pentagon as an example; it is assumed that the centroid of the particle is located at the origin of rectangle coordinates and the x-axis passes through one of its vertices. The size and length of the pentagon are denoted by the notation l. Furthermore, all the vertices P i (i.e., P 1, P 2, P 3, P 4, and P 5 ) of a particle can be symbolized by Equation (2). (2) on the geometric features of the oval particles in Figure 2, the profile of ovals can be simply treated as an extension of ellipses, and their mathematical expression can be described by Equation (1).
where ζ is a tapering coefficient that varies from 0.0 to 1.0, a and b are the semi-axis lengths of the oval in the x-and y-axis directions, respectively, and the ratio of b to a is defined as the aspect ratio of the oval. As shown in Figure 3, a set of regular polygons with different numbers of edges are graphically displayed. According to the study in [18], the geometric topographies of these polygons can be mathematically expressed via the collection of the subsets of vertices and edges. Take the regular pentagon as an example; it is assumed that the centroid of the particle is located at the origin of rectangle coordinates and the x-axis passes through one of its vertices. The size and length of the pentagon are denoted by the notation l. Furthermore, all the vertices Pi (i.e., P1, P2, P3, P4, and P5) of a particle can be symbolized by Equation (2).  Moreover, the five edges of a pentagon are also symbolized by the notation L i (i.e., L 1, L 2, L 3, L 4, and L 5 ), and their detailed expressions can be further mathematically described as: L 1 = {P 1 , P 2 }, L 2 = {P 2 , P 3 }, L 3 = {P 3 , P 4 }, L 4 = {P 4 , P 5 } and L 5 = {P 5 , P 1 }. The relevant geometric information for other kinds of regular polygons has been presented in [18] and is not elaborated again.

Modeling of ITZ around Complex-Shaped Aggregates
In the field of material science, the importance of interphase thickness has been widely acknowledged. For the concrete systems, the thickness of ITZ t ITZ around aggregate particles generally fluctuates in a certain range (e.g., [15 µm, 50 µm] in [31], [50 µm, 200 µm] in [32], [100 µm, 200 µm] in [33], etc.), and its specific value may also be differential at different locations around a single aggregate. Strictly speaking, the values of ITZ thickness in concrete would be closely conditioned by a variety of factors such as the composition, test method, curing condition, mixture ratio, production technology, and so on, which is deemed to be a very complicated issue. However, to accurately describe the quantitative relationship between the ITZ thickness and the structural features (or properties) of materials, the ITZ thickness around the aggregates is generally assumed to be a constant in the numerical modeling of concrete. As the relationship between them is formulized, the corresponding expression can be directly applied by substituting the statistical average of ITZ thickness in real material into it. In this study, the commonly used assumption mentioned above is also adopted for the sake of simplicity. That is to say, all the ITZs in the models are assumed to be the uniform interfacial layers coated on the surface of particles, and the overlap of them can freely occur.
To generate the ITZ with a uniform thickness, various techniques were developed based on the geometric shapes of particles. For the ovals, the outer contour of ITZ around them can be discretized into a limited number of points (i.e., (x ITZ , y ITZ )), and the mathematical expression of these discrete points can be further derived by combining the generic formula of ovals (i.e., Equation (1)) and the unit normal vectors (n x , n y ) of the points on the contour of these particles as follows: ∂F(x, y) ∂y where t ITZ is the thickness of the uniform ITZ around aggregates, (x, y) represent the spatial coordinates of the points on the contour of the oval. For the regular polygons, the outer contour of ITZ around them can be generally subdivided into the part of vertices and the part of edges, as described in [18]. The part of the vertices on the contour of ITZ can be viewed as the circular arcs, which are composed of a limited number of discrete points at a fixed distance from the vertexes. The part of the edges on the contour of ITZ can be treated as the line segments, which are parallel to the original edges of the polygon, and they can be constructed by translating the original edges by a preset distance along the normal line of these edges. of ITZ around them, two common functions (i.e., the EVF function and Fuller function described by Equation (7)) used in the modeling of concrete materials are adopted here to characterize the size gradation (PSG) of these aggregate particles by reference to the literature [34]. Moreover, to uniformly symbolize the sizes of these non-spherical particles with different sizes, the equivalent diameter D eq , which is defined as the diameter of a circle having the same area as that of a 2D complex-shaped particle, is further introduced. Finally, by combining a series of the above preconditions with the RSP algorithm in Figure 4, meso-scale models of polyphase concrete composed of oval fine aggregates, polygonal coarse aggregates, uniform ITZ, and homogeneous cement matrix can be generated in 2D space. The detailed algorithms and technique have been presented in previous studies [18,29,[34][35][36][37]. To avoid duplication, they are not stated in detail in this paper.
where F A (D eq ) is the area-based cumulative probability function, D eq,min and D eq,max are the minimum and maximum equivalent diameters of particles, respectively.  As shown in Figure 5, taking the Fuller function in Equation (7) as an example, one sample of the 2D model of concrete comprising of pentagonal coarse aggregates, oval fine aggregates, uniform ITZ, and cement matrix with the periodic boundary conditions is displayed graphically. In this model, the size length L of a square container is equal to 100 mm. The total fraction ϕagg of aggregate particles is 0.75, and the sand ratio βs (i.e., the area ratio of fine aggregates to the sum of both coarse aggregates and fine aggregates) is 0.40. Besides, all the values of tITZ around the aggregates are set to 100 μm. As shown in Figure 5, taking the Fuller function in Equation (7) as an example, one sample of the 2D model of concrete comprising of pentagonal coarse aggregates, oval fine aggregates, uniform ITZ, and cement matrix with the periodic boundary conditions is displayed graphically. In this model, the size length L of a square container is equal to 100 mm. The total fraction φ agg of aggregate particles is 0.75, and the sand ratio β s (i.e., the  As shown in Figure 5, taking the Fuller function in Equation (7) as an example, one sample of the 2D model of concrete comprising of pentagonal coarse aggregates, oval fine aggregates, uniform ITZ, and cement matrix with the periodic boundary conditions is displayed graphically. In this model, the size length L of a square container is equal to 100 mm. The total fraction ϕagg of aggregate particles is 0.75, and the sand ratio βs (i.e., the area ratio of fine aggregates to the sum of both coarse aggregates and fine aggregates) is 0.40. Besides, all the values of tITZ around the aggregates are set to 100 μm.

Percolation Simulation of ITZ around Polydisperse Aggregates
As of now, a wide variety of percolation models with HCSS networks and numerical techniques (e.g., forest fire models, tree-burning algorithms, etc.) have been developed recently [7,9,15,[17][18][19]26]. In these HCSS models, the criticality of the emergence of the connected paths is generally represented by the critical threshold, which is commonly symbolized by the critical fraction φ agg,c of particulate components. To quantify the ITZ percolation around the mixture of ovals and polygons, a commonly used numerical method in ref. [15] is also adopted in this study, and the detailed procedure for searching the percolation path of ITZ is shown in Figure 6. Firstly, the 2D polyphase model with the periodic boundary condition is generated, as shown in Figure 6a. By taking the left and right edges of the container as an example, the first step is to search for the aggregate particles around which the ITZs can intersect with the left edge and label them with a new tag (i.e., the blue in Figure 6b). Afterward, search further for the aggregate particles around which the ITZs can intersect with the ITZs around the preceding labeled particles and label them with the same tag, as shown in Figure 6c. The above process is performed iteratively until there are no additional particles. Finally, judge whether or not the ITZs around these labeled aggregates can intersect with the right edge of the container. If there exists the interconnectivity among ITZs throughout the whole system, as shown in Figure 6d, the ITZ percolation is deemed to occur in this sample.

Percolation Simulation of ITZ around Polydisperse Aggregates
As of now, a wide variety of percolation models with HCSS networks and numerical techniques (e.g., forest fire models, tree-burning algorithms, etc.) have been developed recently [7,9,15,[17][18][19]26]. In these HCSS models, the criticality of the emergence of the connected paths is generally represented by the critical threshold, which is commonly symbolized by the critical fraction ϕagg,c of particulate components. To quantify the ITZ percolation around the mixture of ovals and polygons, a commonly used numerical method in ref. [15] is also adopted in this study, and the detailed procedure for searching the percolation path of ITZ is shown in Figure 6. Firstly, the 2D polyphase model with the periodic boundary condition is generated, as shown in Figure 6a. By taking the left and right edges of the container as an example, the first step is to search for the aggregate particles around which the ITZs can intersect with the left edge and label them with a new tag (i.e., the blue in Figure 6b). Afterward, search further for the aggregate particles around which the ITZs can intersect with the ITZs around the preceding labeled particles and label them with the same tag, as shown in Figure 6c. The above process is performed iteratively until there are no additional particles. Finally, judge whether or not the ITZs around these labeled aggregates can intersect with the right edge of the container. If there exists the interconnectivity among ITZs throughout the whole system, as shown in Figure  6d, the ITZ percolation is deemed to occur in this sample. In light of the independence of each sample, the occurrence of ITZ percolation in the finite-sized system can be regarded as a Bernoulli probability variable. For a large number of trials, the number of successes could be mathematically represented by the binominal distribution, and the number ratio of the percolating samples Npcl to the total samples Ntotal In light of the independence of each sample, the occurrence of ITZ percolation in the finite-sized system can be regarded as a Bernoulli probability variable. For a large number of trials, the number of successes could be mathematically represented by the binominal distribution, and the number ratio of the percolating samples N pcl to the total samples N total (i.e., P ITZ = N pcl /N total ) is generally defined as the percolation probability P ITZ . A great number of studies [13,15,18,19,37] have clearly shown that the stability of P ITZ is closely dependent on the value of N total . To improve the statistical stability of P ITZ , N total is set to be 2000 here by referring to the description in [18]. Furthermore, as the probability P ITZ for the finite-sized systems with different L and φ agg is obtained, the critical percolation threshold φ agg,c can be further determined by the curve-fitting method with the following formula (i.e., Equation (8)).
where φ agg, is the relative content of aggregate particles, φ agg,c is the critical percolation threshold of ITZ, and ∆(L) is the percolation transition width. According to the study by Pan et al. [26], the size length L of percolation models would mainly affect the value of ∆(L), but has no great influence on the critical percolation threshold φ agg,c . By reference to the findings (i.e., L ≥ 5D eq,max ) in [19] and the criterion for determining the representative volume element (REV) in [38] (i.e., the RVE size should be at least 3~5 times of the maximum particle diameter), the size lengths of all the generated models here are set to not less than 5 times of the maximum size D eq,max,P of polygonal coarse aggregates (i.e., L/D eq,max,P ≥ 5.0).

Results and Discussion
In the section, a series of parametric analyses about the effect of aggregate features on the critical threshold φ agg,c of ITZ percolation are conducted. Based on the simulated results of φ agg,c , an approximately analytical function is further developed to fast evaluate the ITZ percolation in the 2D case of polyphase composite systems.

Effect of Particle Shapes on the ITZ Percolation
By assuming the coarse aggregates as regular pentagons, the sensitivity of the critical threshold φ agg,c to the shape of oval fine aggregates is studied first in Figure 7. The curved surfaces clearly show that the threshold φ agg,c is an apparently ζ-and b/a-dependent argument, which possesses the obviously different trends with the parameters ζ and b/a. On the one hand, as ζ increases from 0.0 to 1.0, φ agg,c basically shows a slightly decreasing trend under a constant b/a. On the other hand, when the coefficient ζ is fixed, all the φ agg,cb/a curves could be subdivided into two stages (i.e., the upward stage and the downward stage) with the increase of b/a in [0.2, 2.5], and the inflection points are basically located at b/a = 1.0. In the upward stage, the value of φ agg,c increases rapidly with the increasing b/a from 0.2 to 1.0 in the parabolic manner. However, with the further growth of b/a, φ agg,c presents a very slowly descending trend. Overall, the maximal threshold φ agg,c appears at ζ = 0.0 and b/a = 1.0, which indicates that the ITZ percolation paths are least likely to be formed if the fine aggregates are assumed to be circular. In [29], the variation curves of the percolation threshold ε c,mm of porous systems composed of oval pores with b/a in [0.2, 5.0] and ζ in [0.0, 0.5] were presented. By comparison, one can clearly see that the variation of φ agg,c with b/a and ζ here is commendably consistent with these results. These comparisons indirectly indicate the reliability of the variation of φ agg,c .
In the literature [18], the aggregate shapes are uniformly symbolized by the circularity c 2D (i.e., the perimeter ratio of circles to non-circular particles having the same area), and the critical threshold φ agg,c is formulized as a function of c 2D . By referring to this strategy, the circularities c 2D,O of ovals with b/a in [0.2, 2.5] and ζ in [0.0, 1.0] are calculated, and the sensitivity of φ agg,c to the circularity c 2D,O is further graphically presented based on the results of φ agg,c in Figure 7. It can be seen that all the data points in Figure 8 are roughly converging to a quadratic increasing curve, which means that the smaller the circularity c 2D,O of fine aggregates is, the more easily the percolation paths of ITZ around them are formed. The above trend indicates that the circularity can be well used to represent the quantitative relationship between the ITZ percolation threshold and the oval shapes. This provides the foundation for the derivation of the analytical model of φ agg,c in the following sections. ment, which possesses the obviously different trends with the parameters ζ and b/a. On the one hand, as ζ increases from 0.0 to 1.0, ϕagg,c basically shows a slightly decreasing trend under a constant b/a. On the other hand, when the coefficient ζ is fixed, all the ϕagg,cb/a curves could be subdivided into two stages (i.e., the upward stage and the downward stage) with the increase of b/a in [0.2, 2.5], and the inflection points are basically located at b/a = 1.0. In the upward stage, the value of ϕagg,c increases rapidly with the increasing b/a from 0.2 to 1.0 in the parabolic manner. However, with the further growth of b/a, ϕagg,c presents a very slowly descending trend. Overall, the maximal threshold ϕagg,c appears at ζ = 0.0 and b/a = 1.0, which indicates that the ITZ percolation paths are least likely to be formed if the fine aggregates are assumed to be circular. In [29], the variation curves of the percolation threshold εc,mm of porous systems composed of oval pores with b/a in [0.2, 5.0] and ζ in [0.0, 0.5] were presented. By comparison, one can clearly see that the variation of ϕagg,c with b/a and ζ here is commendably consistent with these results. These comparisons indirectly indicate the reliability of the variation of ϕagg,c. In the literature [18], the aggregate shapes are uniformly symbolized by the circularity c2D (i.e., the perimeter ratio of circles to non-circular particles having the same area), and the critical threshold ϕagg,c is formulized as a function of c2D. By referring to this strategy, the circularities c2D,O of ovals with b/a in [0.2, 2.5] and ζ in [0.0, 1.0] are calculated, and the sensitivity of ϕagg,c to the circularity c2D,O is further graphically presented based on the results of ϕagg,c in Figure 7. It can be seen that all the data points in Figure 8 are roughly converging to a quadratic increasing curve, which means that the smaller the circularity c2D,O of fine aggregates is, the more easily the percolation paths of ITZ around them are  In Figure 9, the sensitivity of ITZ percolation to the shape of polygonal coarse aggregates is further studied. The results clearly show that no matter what the shape of fine aggregates is, ϕagg,c always possesses a similar trend with the shape shifting of coarse aggregates (i.e., Circle > Decagon > Octagon > Hexagon > Pentagon > Square > Triangle). By combining with the circularity c2D,P of these polygons, the value of ϕagg,c here also shows the downward trend with the increasing c2D,P and the attenuation amplitude of ϕagg,c is largely dependent on the relative difference of c2D,P for different polygons. However, by contrast with the dataset in [18], it can be found that the impact degree of aggregate shape (i.e, c2D,O or c2D,P) on the value of ϕagg,c here is significantly weakened, which indicates that the analytical formula in the previous paper [18] may not be directly applicable. Overall, the ITZ percolation would be influenced by the synergy of both fine aggregate shapes and coarse aggregate shapes. It is very necessary to distinguish the geometric difference between fine aggregate and coarse aggregate in the modeling of concrete. In Figure 9, the sensitivity of ITZ percolation to the shape of polygonal coarse aggregates is further studied. The results clearly show that no matter what the shape of fine aggregates is, φ agg,c always possesses a similar trend with the shape shifting of coarse aggregates (i.e., Circle > Decagon > Octagon > Hexagon > Pentagon > Square > Triangle). By combining with the circularity c 2D,P of these polygons, the value of φ agg,c here also shows the downward trend with the increasing c 2D,P and the attenuation amplitude of φ agg,c is largely dependent on the relative difference of c 2D,P for different polygons. However, by contrast with the dataset in [18], it can be found that the impact degree of aggregate shape (i.e, c 2D,O or c 2D,P ) on the value of φ agg,c here is significantly weakened, which indicates that the analytical formula in the previous paper [18] may not be directly applicable. Overall, the ITZ percolation would be influenced by the synergy of both fine aggregate shapes and coarse aggregate shapes. It is very necessary to distinguish the geometric difference between fine aggregate and coarse aggregate in the modeling of concrete. largely dependent on the relative difference of c2D,P for different polygons. However, by contrast with the dataset in [18], it can be found that the impact degree of aggregate shape (i.e, c2D,O or c2D,P) on the value of ϕagg,c here is significantly weakened, which indicates that the analytical formula in the previous paper [18] may not be directly applicable. Overall, the ITZ percolation would be influenced by the synergy of both fine aggregate shapes and coarse aggregate shapes. It is very necessary to distinguish the geometric difference between fine aggregate and coarse aggregate in the modeling of concrete. Figure 9. Effect of the geometric shapes of polygonal particles on the ITZ percolation. Figure 9. Effect of the geometric shapes of polygonal particles on the ITZ percolation.

Effect of Particle Sizes on the ITZ Percolation
Size-polydispersity of particles inevitably exists in most granular materials, and a number of studies about the ITZ percolation involving the polysized aggregates have been conducted recently [13,15,18,19,26], but the size distribution difference between fine aggregates and coarse aggregates is often overlooked. In the section, the dependence of critical threshold φ agg,c on the size distributions of both fine and coarse aggregates is separately analyzed. Figure 10a,b show the sensitivity of φ agg,c to the maximum diameter D eq,max,O of fine aggregates and maximum diameter D eq,max,P of coarse aggregates, respectively. It can be seen that all the curves of φ agg,c -t ITZ possess a remarkable inverted parabolic downward trend with the increasing t ITZ in [0.03 mm, 0.11 mm]. This is mainly due to the fact that the thicker the ITZ thickness is, the greater the collision probability of ITZ around different aggregates is, which in turn promotes the formation of ITZ percolation paths in concrete. Besides, when the value of D eq,min,O or D eq,min,P is held constant, the threshold φ agg,c shows the various degrees of growth with the increase of D eq,max,O or D eq,max,P . That is to say, the smaller the value of D eq,max,O or D eq,max,P is, the more easily the ITZ percolation paths are formed. By comparison, it is not difficult to find that the dependence of φ agg,c on the size D eq,max,P is significantly lower than the size D eq,max,O . Hence, for the sake of simplicity, the effect of the size ranges of coarse aggregates on the ITZ percolation can be negligible in some cases.

Effect of Particle Sizes on the ITZ Percolation
Size-polydispersity of particles inevitably exists in most granular materials, and a number of studies about the ITZ percolation involving the polysized aggregates have been conducted recently [13,15,18,19,26], but the size distribution difference between fine aggregates and coarse aggregates is often overlooked. In the section, the dependence of critical threshold ϕagg,c on the size distributions of both fine and coarse aggregates is separately analyzed. Figure 10a,b show the sensitivity of ϕagg,c to the maximum diameter Deq,max,O of fine aggregates and maximum diameter Deq,max,P of coarse aggregates, respectively. It can be seen that all the curves of ϕagg,c-tITZ possess a remarkable inverted parabolic downward trend with the increasing tITZ in [0.03 mm, 0.11 mm]. This is mainly due to the fact that the thicker the ITZ thickness is, the greater the collision probability of ITZ around different aggregates is, which in turn promotes the formation of ITZ percolation paths in concrete. Besides, when the value of Deq,min,O or Deq,min,P is held constant, the threshold ϕagg,c shows the various degrees of growth with the increase of Deq,max,O or Deq,max,P. That is to say, the smaller the value of Deq,max,O or Deq,max,P is, the more easily the ITZ percolation paths are formed. By comparison, it is not difficult to find that the dependence of ϕagg,c on the size Deq,max,P is significantly lower than the size Deq,max,O. Hence, for the sake of simplicity, the effect of the size ranges of coarse aggregates on the ITZ percolation can be negligible in some cases. By using Equation (7), the effect of PSG of both fine aggregates and coarse aggregates on the ITZ percolation is further shown in Figure 10c. One can clearly see that the PSG of fine aggregates has a significant influence on the critical threshold ϕagg,c. The value of ϕagg,c for EVF gradation is much smaller than the corresponding threshold for Fuller gradation. This is largely due to the fact that the number of small-sized aggregates in the system for By using Equation (7), the effect of PSG of both fine aggregates and coarse aggregates on the ITZ percolation is further shown in Figure 10c. One can clearly see that the PSG of fine aggregates has a significant influence on the critical threshold φ agg,c . The value of φ agg,c for EVF gradation is much smaller than the corresponding threshold for Fuller gradation. This is largely due to the fact that the number of small-sized aggregates in the system for EVF gradation is generally larger than that for Fuller gradation. The above trend of φ agg,c here is in good agreement with the existing findings for the mono-shaped particle system in [15,18,19,26]. In contrast to that, the PSG of coarse aggregates has little effect on the ITZ percolation around them. This phenomenon has not yet been reported in previous literature and needs to be given enough consideration.
According to the study in [18], the influence of all the size-polydispersities of aggregate particles on the ITZ percolation φ agg,c can be attributed to the aggregate fineness, which is commonly symbolized by the specific surface area C A . Based on this view, the average specific surface area C avg A of both oval aggregates and polygonal aggregates is derived here to represent the fineness of aggregate particles in concrete, as expressed by Equations (9)- (11).
C A,P = D eq,max,P D eq,min,P C P D eq,P · f N D eq,P dD eq,P D eq,max,P D eq,min,P A P D eq,P · f N D eq,P dD eq,P , (11) where C A,O and C A,P are the specific surface areas of oval fine aggregates and polygonal coarse aggregates, C O (D eq,O ), A O (D eq,O ) and f N (D eq,O ) are the perimeter, area and numberbased probability function for oval fine aggregates, respectively; C P (D eq,P ), A P (D eq,P ), and f N (D eq,P ) are the perimeter, area, and number-based probability function for polygonal coarse aggregates, respectively. For the Fuller function and EVF function in Equation (7), the corresponding numberbased probability functions can be expressed as below: By substituting Equation (12) into Equations (10) and (11), the theoretical solution of C A for different types of particles can be further derived, and the detailed formulas (i.e., Equation (13)) are expressed as a function of c 2D .
Finally, by incorporating Equation (13) into Equation (9), the average specific surface area C avg A of aggregate particles for all the cases in Figure 10 is obtained, and the comparisons of them are shown in Figure 11a-c. From these figures, the following trends can be found: (a) as shown in Figure 11a,b, with the increase of D eq,max,O or D eq,max,P , the value of C avg A is also increased significantly, and its sensitivity to the size range of fine aggregates is much bigger than that for coarse aggregates; (b) The PSG of coarse aggregates has little influence on the value of C avg A , however, when the PSG of coarse aggregates is fixed, C avg A for the fine aggregates with EVF gradation is significantly greater than the value of C avg A with Fuller gradation, as shown in Figure 11c. By combining with the curves in Figure 10, a relatively consistent trend can be provided (i.e., the greater the average specific surface area of aggregate particles is, the smaller the percolation threshold of ITZ φ agg,c tends to be), which also provides a direction for the derivation of the model of φ agg,c .

Effect of Sand Ratio on the ITZ Percolation
The sand ratio βs is an important factor that affects the workability and mechanical properties of concrete materials. Under the assumption of Deq,O = Fuller 0.15-4.75 mm and Deq,P = Fuller 5-16 mm, the sensitivity of the critical threshold ϕagg,c to the sand ratio βs is investigated in Figure 12. One can see that all the curves of ϕagg,c-βs show an apparently linear decreasing tendency, which is dependent on the specific thickness of ITZ. That is to say, the greater the proportion of fine aggregates in concrete is, the more easily the ITZ percolation path is formed. Actually, the effect of βs on ϕagg,c can also be mapped to the variation of aggregate fineness. As the sand ratio βs increases, the number of fine aggregates is also on the rise, which in turn leads to the rapid growth of the average specific surface area of aggregate particles, as shown in Figure 13.

Effect of Sand Ratio on the ITZ Percolation
The sand ratio β s is an important factor that affects the workability and mechanical properties of concrete materials. Under the assumption of D eq,O = Fuller 0.15-4.75 mm and D eq,P = Fuller 5-16 mm, the sensitivity of the critical threshold φ agg,c to the sand ratio β s is investigated in Figure 12. One can see that all the curves of φ agg,c -β s show an apparently linear decreasing tendency, which is dependent on the specific thickness of ITZ. That is to say, the greater the proportion of fine aggregates in concrete is, the more easily the ITZ percolation path is formed. Actually, the effect of β s on φ agg,c can also be mapped to the variation of aggregate fineness. As the sand ratio β s increases, the number of fine aggregates is also on the rise, which in turn leads to the rapid growth of the average specific surface area of aggregate particles, as shown in Figure 13.

Effect of Sand Ratio on the ITZ Percolation
The sand ratio βs is an important factor that affects the workability and mechanical properties of concrete materials. Under the assumption of Deq,O = Fuller 0.15-4.75 mm and Deq,P = Fuller 5-16 mm, the sensitivity of the critical threshold ϕagg,c to the sand ratio βs is investigated in Figure 12. One can see that all the curves of ϕagg,c-βs show an apparently linear decreasing tendency, which is dependent on the specific thickness of ITZ. That is to say, the greater the proportion of fine aggregates in concrete is, the more easily the ITZ percolation path is formed. Actually, the effect of βs on ϕagg,c can also be mapped to the variation of aggregate fineness. As the sand ratio βs increases, the number of fine aggregates is also on the rise, which in turn leads to the rapid growth of the average specific surface area of aggregate particles, as shown in Figure 13.

Formulation and Verification
Through a series of in-depth analyses and comparisons, it can be found that the ITZ thickness, aggregate fineness, and aggregate shape can be deemed to be the main contributing parameters that affect the ITZ percolation when the synergetic effects of aggregate shapes and aggregate sizes are taken into consideration. To quantitatively characterize the sensitivity of the critical percolation threshold ϕagg,c to all of the aforementioned factors, a new argument η (i.e., Equation (14)) is proposed here by referring to the form of the proposed parameter (i.e., −ln(tITZCA)) in [18] and all the simulated results of ϕagg,c in Figures 7, 9, 10 and 12 are plotted as a function of η, as shown in Figure 14. The curve of ϕagg,c-η in Figure 14 clearly shows that ϕagg,c possesses an apparently monotonical uptrend with the increase of η in a quadratic manner. By the linear regression

Formulation and Verification
Through a series of in-depth analyses and comparisons, it can be found that the ITZ thickness, aggregate fineness, and aggregate shape can be deemed to be the main contributing parameters that affect the ITZ percolation when the synergetic effects of aggregate shapes and aggregate sizes are taken into consideration. To quantitatively characterize the sensitivity of the critical percolation threshold φ agg,c to all of the aforementioned factors, a new argument η (i.e., Equation (14)) is proposed here by referring to the form of the proposed parameter (i.e., −ln(t ITZ C A )) in [18] and all the simulated results of φ agg,c in Figures 7, 9, 10 and 12 are plotted as a function of η, as shown in Figure 14.
Materials 2023, 16, x FOR PEER REVIEW 14 of 19 Figure 13. Effect of the sand ratio βs on the average specific surface area of particles for the cases in Figure 12.

Formulation and Verification
Through a series of in-depth analyses and comparisons, it can be found that the ITZ thickness, aggregate fineness, and aggregate shape can be deemed to be the main contributing parameters that affect the ITZ percolation when the synergetic effects of aggregate shapes and aggregate sizes are taken into consideration. To quantitatively characterize the sensitivity of the critical percolation threshold ϕagg,c to all of the aforementioned factors, a new argument η (i.e., Equation (14)) is proposed here by referring to the form of the proposed parameter (i.e., −ln(tITZCA)) in [18] and all the simulated results of ϕagg,c in Figures 7, 9, 10 and 12 are plotted as a function of η, as shown in Figure 14. The curve of ϕagg,c-η in Figure 14 clearly shows that ϕagg,c possesses an apparently monotonical uptrend with the increase of η in a quadratic manner. By the linear regression The curve of φ agg,c -η in Figure 14 clearly shows that φ agg,c possesses an apparently monotonical uptrend with the increase of η in a quadratic manner. By the linear regression analysis method, a numerical fitting formula is obtained in Figure 14, and the approximately empirical expression of φ agg,c (i.e., Equation (15)) can be further derived by incorporating the argument η into this quadratic formula. To validate the reliability of the proposed function of φ agg,c , the comparison of the simulated values of φ agg,c in some extra cases with the predicted approximation obtained by Equation (15) is shown in Figure 15a. It can be seen that all the curves of φ agg,c -t ITZ here are in excellent agreement with the simulated results. Figure 15b further displays the quantitative comparison of the theoretical solution of φ agg,c with the numerical results reported by Zheng et al. [15], and the approximate agreement between them can also be observed clearly. All these comparisons suggest that the proposed analytical formula of φ agg,c can reproduce the results of φ agg,c with the satisfying accuracy, and it can be well used to quickly evaluate the ITZ percolation in 2D polyphase systems. analysis method, a numerical fitting formula is obtained in Figure 14, and the approximately empirical expression of ϕagg,c (i.e., Equation (15)) can be further derived by incorporating the argument η into this quadratic formula. To validate the reliability of the proposed function of ϕagg,c, the comparison of the simulated values of ϕagg,c in some extra cases with the predicted approximation obtained by Equation (15) is shown in Figure 15a. It can be seen that all the curves of ϕagg,c-tITZ here are in excellent agreement with the simulated results. Figure 15b further displays the quantitative comparison of the theoretical solution of ϕagg,c with the numerical results reported by Zheng et al. [15], and the approximate agreement between them can also be observed clearly. All these comparisons suggest that the proposed analytical formula of ϕagg,c can reproduce the results of ϕagg,c with the satisfying accuracy, and it can be well used to quickly evaluate the ITZ percolation in 2D polyphase systems.
As shown in Figure 16, it is assumed that the ITZ thickness tITZ around aggregates is fixed to be 0.05 mm. The size distributions of aggregate particles are respectively set to Deq = Fuller 0.15-10 mm and Deq = Fuller 0.15-16 mm. As we know from the calculation, the corresponding sand ratio βs in them can be separately equal to 0.50 and 0.65. Take the above cases as the study object, and the predicted approximation of ϕagg,c by the proposed functions (i.e., Equations (9), (13) and (15)) here is presented, and compared with the analytical results predicted by the formula in [18]. It can be seen that with the decrease of c2D,O in [0.5, 1.0] or c2D,P in [0.75, 1.0], the critical threshold of ITZ ϕagg,c also decreases apparently. Moreover, when the circularities of both fine aggregates and coarse aggregates are identical (i.e., c2D,O = c2D,P), the proposed functions generally provide slightly larger results of ϕagg,c. than the results in [18], and the smaller the particle circularity is, the larger the relative error between them would be. This can be blamed on two main reasons. The first one is that the analytical formula in [18] was mainly derived by analyzing the large bodies of data for the ITZ percolation around circles with its size less than 5 mm and the small amount of data for the models of polygons. Literally, the published expression may be more applicable for the 2D mortar systems that composed the circular aggregates. Once it is extended to the concrete systems, including the non-circular aggregates with sizes larger than 5 mm, its effectiveness and accuracy may be limited to some extent. The second reason is that, compared with the existing literature, the size gradations of both coarse aggregate and fine aggregates in this work are separately expressed, which is more practical than the single function. Besides, the size ranges of aggregate particles used in this work are also much larger and wider than before. Simultaneously, some other important As shown in Figure 16, it is assumed that the ITZ thickness t ITZ around aggregates is fixed to be 0.05 mm. The size distributions of aggregate particles are respectively set to D eq = Fuller 0.15-10 mm and D eq = Fuller 0.15-16 mm. As we know from the calculation, the corresponding sand ratio β s in them can be separately equal to 0.50 and 0.65. Take the above cases as the study object, and the predicted approximation of φ agg,c by the proposed functions (i.e., Equations (9), (13) and (15)) here is presented, and compared with the analytical results predicted by the formula in [18]. It can be seen that with the decrease of c 2D,O in [0.5, 1.0] or c 2D,P in [0.75, 1.0], the critical threshold of ITZ φ agg,c also decreases apparently. Moreover, when the circularities of both fine aggregates and coarse aggregates are identical (i.e., c 2D,O = c 2D,P ), the proposed functions generally provide slightly larger results of φ agg,c. than the results in [18], and the smaller the particle circularity is, the larger the relative error between them would be. This can be blamed on two main reasons. The first one is that the analytical formula in [18] was mainly derived by analyzing the large bodies of data for the ITZ percolation around circles with its size less than 5 mm and the small amount of data for the models of polygons. Literally, the published expression may be more applicable for the 2D mortar systems that composed the circular aggregates. Once it is extended to the concrete systems, including the non-circular aggregates with sizes larger than 5 mm, its effectiveness and accuracy may be limited to some extent. The second reason is that, compared with the existing literature, the size gradations of both coarse aggregate and fine aggregates in this work are separately expressed, which is more practical than the single function. Besides, the size ranges of aggregate particles used in this work are also much larger and wider than before. Simultaneously, some other important factors (e.g., the aggregate shapes and sand ratio) are also considered in this work. All of the above conditions make the predicted results here more accurate and comprehensive. Through the above figures, it can also be seen that the proposed model here is applicable not only to the granular systems of mono-shaped particles (c 2D,O = c 2D,P ), but also to the systems involving the effect of the shape diversity (c 2D,O = c 2D,P )and relative content (β s ) of different particles. Overall, the model here is bound to have a much broader range of applications in the 2D case compared with the published formula from the previous paper. factors (e.g., the aggregate shapes and sand ratio) are also considered in this work. All of the above conditions make the predicted results here more accurate and comprehensive. Through the above figures, it can also be seen that the proposed model here is applicable not only to the granular systems of mono-shaped particles (c2D,O = c2D,P), but also to the systems involving the effect of the shape diversity (c2D,O ≠ c2D,P)and relative content (βs) of different particles. Overall, the model here is bound to have a much broader range of applications in the 2D case compared with the published formula from the previous paper. shapes (i.e., c2D,P and c2D,O) can be deemed to the contributing factors for the percolation properties of ITZ around the 2D poly-shaped poly-sized aggregates. However, given that all the work in this paper is conducted based on the 2D simplified granular systems, the proposed formulas here can only be applicable to the 2D granular structures. That is to say, it is still difficult to be directly applied in engineering. But, the variation trends of the ITZ percolation with the influencing factors should be basically consistent in both the 2D and 3D cases. Hence, the work in this paper can still play a positive role in the development of similar research in 3D composite systems. Finally, strictly speaking, the geometric features of aggregate particles in actual concrete materials are much more complex and diverse than that in the models here. By reference to the proposed formula (i.e., Equation (15)) here, an analytical formula of ITZ percolation for the more complex systems is preliminarily suggested, as expressed by Equation (16). The specific values of the parameters (i.e., a1, a2, and a3) and the feasibility of Equation (16) need to be further studied. It is hoped that this hypothesis may provide some help for others in their research.
where a1, a2 and a3 are three coefficients, CA,i and c2D,i are the specific surface area and circularity of the ith types of particles, respectively, βi is defined as the relative fraction of the ith types of particles and the sum of β1 + β2 + … βi … + βN-1 + βN is exactly equal to 1.0, in which N is the total number of particle species in the system. According to the curves of φ agg,c in Figures 7-14, it can be found that all the aforementioned factors including ITZ thickness (t ITZ ), aggregate fineness (C avg A ) and aggregate shapes (i.e., c 2D,P and c 2D,O ) can be deemed to the contributing factors for the percolation properties of ITZ around the 2D poly-shaped poly-sized aggregates. However, given that all the work in this paper is conducted based on the 2D simplified granular systems, the proposed formulas here can only be applicable to the 2D granular structures. That is to say, it is still difficult to be directly applied in engineering. But, the variation trends of the ITZ percolation with the influencing factors should be basically consistent in both the 2D and 3D cases. Hence, the work in this paper can still play a positive role in the development of similar research in 3D composite systems.
Finally, strictly speaking, the geometric features of aggregate particles in actual concrete materials are much more complex and diverse than that in the models here. By reference to the proposed formula (i.e., Equation (15)) here, an analytical formula of ITZ percolation for the more complex systems is preliminarily suggested, as expressed by Equation (16). The specific values of the parameters (i.e., a 1 , a 2 , and a 3 ) and the feasibility of Equation (16) need to be further studied. It is hoped that this hypothesis may provide some help for others in their research.
where a 1 , a 2 and a 3 are three coefficients, C A,i and c 2D,i are the specific surface area and circularity of the ith types of particles, respectively, β i is defined as the relative fraction of the ith types of particles and the sum of β 1 + β 2 + . . . β i . . . + β N−1 + β N is exactly equal to 1.0, in which N is the total number of particle species in the system.

Conclusions
In this paper, the sensitivity of ITZ percolation to the features of both coarse aggregates and fine aggregates is studied in a 2D case. To highlight the effect of the shape and size difference between coarse aggregates and fine aggregates in a more realistic way, they are respectively simplified as polygons and ovals here based on the real morphologies of these aggregates. Moreover, the size gradations of them are also expressed by two separate functions instead of the single contiguous function in the literature. By coupling the generated 2D polyphase models of concrete systems with the continuum percolation, the percolation behavior of the uniform ITZ around aggregates is simulated and the effect of the polydispersity of aggregate shapes and sizes on the critical threshold φ agg,c is quantitatively analyzed, respectively. Numerical variation of the threshold φ agg,c reveals that: (a) the ITZ thickness (t ITZ ), average aggregate fineness (C avg A ), coarse aggregate shape (c 2D,P ), and fine aggregate shape (c 2D,O ) are four important factors to the ITZ percolation, (b) the sand ratio β s is also an important parameter to the ITZ percolation and its effect can be mapped to the aggregate fineness (i.e., the larger the ratio β s is, the larger the average specific surface of aggregates would be, which in turn leads to the smaller values of φ agg,c ).
By adopting the linear regression analysis, an approximately analytical expression of φ agg,c (i.e., Equations (9), (13) and (15)) is further proposed and its applicability is verified by comparing with the numerical results, which indicates that the proposed formulas here can reproduce the results of φ agg,c with the satisfying accuracy. Moreover, by comparing the proposed model with the existing prediction expression, it can be found that the model can produce slightly larger results for the ITZ percolation threshold φ agg,c , which may be closer to reality. Finally, it can not only be applied to the mortar or concrete but also to some other similar granular systems in the 2D case.  x-axis component of unit normal vector D eq,min,O minimum equivalent diameter of oval n y y-axis component of unit normal vector D eq,max,P maximum equivalent diameter of polygon β s sand ratio D eq,min,P minimum equivalent diameter of polygon F A area-based cumulative probability function C A specific surface area of particle D eq equivalent diameter of particle C A,O, C A,P specific surface area of ovals and polygons, respectively D eq,max maximum equivalent diameter of particle C A,i specific surface area of the ith type particle D eq,min minimum equivalent diameter of particle C avg A average surface area of both fine and coarse aggregate particles D eq,O equivalent diameter of oval C O , A O perimeter and area of oval, respectively D eq,P equivalent diameter of polygon C P , A P perimeter and area of polygon, respectively N pcl number of the percolating samples f N number-based probability function of particles N total number of the total samples β i fraction of the ith types of particles P ITZ percolation probability a 1 , a 2 , a 3 three coefficients L size length of 2D square models