A 3D multi-phase meso-scale model for modelling coupling of damage and transport properties of concrete

A 3D multi-phase meso-scale model is proposed to analyse the evolution of damage and transport properties in concrete. This model includes five phases: aggregates, mortar, voids, cracking and interfacial transition zones (ITZs) between aggregates and mortar. They are generated using measurable characteristics, such as shape, size and volume fraction of aggregates and voids. Mechanical part of the model is validated through comparing the predicted and experimental stress-strain curves and cracking patterns. Transport part of the model is validated by the comparison of predicted and experimental diffusivities of chloride ions in concrete. A new interface element type is also proposed to represent cracking in transport analysis which can simulate transport not only in normal directions but also in tangential directions. The effects of parameters relevant to aggregates, voids, ITZs and cracking patterns on the mechanical and transport responses of concrete are also analysed. The derived quali- tative understanding of the relation between meso-structural features and macroscopic properties will help future advances in modelling and formulation of experimental programmes to improve the calibration and validation of meso-scale models. The proposed sequential mechanical-transport coupling model is applicable to a large class of composite materials with appropriate selection of phases and their properties. mechanical of ITZs from


Introduction
During service, concrete is interacting with environment and detrimental chemical species can penetrate the bulk, leading to degradation of properties via chemical reactions and consequently reduced durability. The transport of detrimental species, via diffusion and/or convection, will be accelerated when mechanical damage and cracking develops. Understanding the mechanical effects on transport in concrete appears to be critical for performance assessment throughout concrete service life and durability evaluation.
Concrete is a complex composite with sub-continuum structures at multiple length scales. The largest length-scale with observable heterogeneity is traditionally referred to as the concrete meso-scale. The structure at this scale is a multiphase material consisting of voids, aggregates, mortar matrix, and interfacial transition zone (ITZ). Being closest to the component scale, the meso-structure plays a critical role in the observable macroscopic behaviour of concrete [1]. Therefore, the experiments and models established in this mesoscale have drawn more and more attention in recent years.
In the multiphase composition of concrete, ITZ is very thin layers of material surrounding the aggregates [2]. Critical feature of ITZ is that on average they are more porous than mortar, with porosity gradually decreasing from a maximum at the aggregate surface to a minimum reaching the mortar. This rather diffuse definition makes the determination of ITZs thickness approximate, but it is typically accepted to be in the range of 10-80 μm [3,4]. Because of their higher porosity, ITZs have stronger impact on the concrete durability than the other solid phases. Higher porosity translates into weaker mechanical response and larger transport coefficients, i.e. permeability or diffusivity. It makes the chloride ions diffusion coefficients from 2 to 12 times higher than the matrix, and reduction of 0%-50% both the compressive and tensile strength of those the mortar [5][6][7]. The weaker mechanical response of ITZs renders ITZ regions of least resistance to damage/crack initiation, while the larger transport coefficients of ITZs render ITZ regions of least resistance to ingress of detrimental species, gaseous or liquid, from the environment. However, a number of previous works did not consider ITZs as a separate phase and assumed perfect bonding between aggregates and mortar [8,9]. This appears to be an oversimplification due to different mechanical and transport behaviour of ITZs from the mortar. These differences can be accounted for by modelling ITZs with either zero-thickness cohesive interface elements [1,7], or with continuum elements of prescribed thickness [2,10]. When accurate calibration of parameters is made, both approaches have been shown to work well in predicting macroscopic stress strain behaviour of concrete, with negligible differences in the pre-peak and somewhat larger differences in the post-peak behaviour. However, a consideration of enhanced transport of detrimental species along ITZ requires the adoption of the second approachcontinuum elements with prescribed thickness. An exact representation of the real ITZ thickness with solid elements at the meso-scale would lead to too fine meshes and prohibitive analysis times at present. Therefore, in previous numerical studies [2,6,10,11] the ITZ thickness was assumed to be between 0.1 mm and 1 mm. Kim and Abu Al-Rub [10] found that the change of ITZ thickness from 0.1 mm to 0.8 mm had little effect on both the pre-peak and post-peak behaviour of concrete. Based on these past results, the present study uses 0.1 mm thick solid elements to represent ITZ as a separate phase with tensile and compressive strength approximately 75% of the mortar's as suggested in Lu and Tu [12].
A number of meso-scale models, considering the meso-structure and the local failure mechanisms in concrete, have been proposed in the last few decades, aiming at better understanding the effect of structures on the longer-scale response of concrete [13,14]. Continuum based finite element models are the main approaches employed in the literature [15,16], mostly based on the damage-plasticity model [10,17] and the cohesive zone model [18][19][20][21]. The damage-plasticity model provides a general capability for modelling concrete and other quasi-brittle materials by assuming damaged elasticity in combination with tensile and compressive plasticity. In the cohesive zone models, cohesive interface elements are pre-or dynamically inserted into the finite element mesh so that realistic crack patterns can be simulated [7]. Alternative approaches have also been developed, such as discrete element models [22] and lattice models [23]. One challenge with these is the calculation of equivalent model parameters, i.e. force-displacement relations need to be calibrated by macroscopic stress-strain data and crack diffusion coefficient needs to be obtained by trial-and-error approaches. This is avoided when using continuum based FE model. The nonlinear finite element simulations conducted in this work for damage and failure analysis of concrete are based on the damage-plasticity model [24]. Aggregates are assumed to be elastic, and different softening constitutive laws for mortar and ITZs are specified. This allows damage to initiate and localise in crack-like regions through the mortar and along the ITZs.
The generation of actual 3D model with distributed random aggregates is not a trivial task, and most of the meso-scale modelling of concrete has been conducted with 2D samples. Further, the majority of past works treat concrete as a two phase (mortar, aggregates) or three phase (mortar, aggregates and interfaces) material at the meso-scale [25][26][27]. However, X-ray Computed Tomography images of concrete show clearly that voids exist in the concrete at this scale [28]. Therefore, volume fraction, size, gradation, spatial distribution of aggregates and voids should be considered as important for adequate meso-structure description. Several methods have been used for generation of 3D concrete meso-structures: Voronoi tessellation [29], lattice model [30], and Monte Carlo random sampling [31], etc. In these models, however, volume fraction, size distribution and shape of aggregates are not explicitly considered. In this paper, 3D heterogeneous concrete specimens with randomly distributed spherical, ellipsoidal, polyhedral aggregates and spherical voids are generated by a three-level hierarchical method [7,32], where aggregates and voids obey prescribed distributions of size, shape and volume fraction.
In the past decades, a number of numerical methods such as finite element, finite difference, random walk and Lattice Boltzmann for investigating transport behaviour or mechanical response in sound and cracked cement-based materials have already been developed [32][33][34][35][36]. However, little numerical models have been developed to couple the 3D mechanical damage and transport process in concrete which simulate the aggregates, ITZ, mortar matrix and voids explicitly to the authors' knowledge. Microstructure characterisations of cement and concrete show a wide range of micro-and millimetre-pore systems in concrete. The actual flow mechanisms are considerably more complex than what can be simulated by conventional dual-porosity models and Darcy's law. A model capturing multiple porosity and corresponding transport properties for different concrete phases is required to improve the understanding of complex flow mechanism in concrete. In this work the concrete pore space is divided into five porosity systems: aggregates, voids, cracks, ITZs and mortar matrix. Aggregates are impermeable for diffusion, whereas other phases control the diffusion process. This new built multiple-porosity model allows for improved representation of the complex flow mechanisms.
The aim of this study is to clarify the effects of concrete mesostructure characteristics on its macroscopic behaviour including a practical approach for coupling the evolutions of mechanical damage and transport coefficients. It is anticipated that the work will inspire further developments and contribute to the design of improved durability concrete. Effective and efficient simulation of 3D damage evolution and failure of concrete at the meso-scale is now possible, due to the recent advances in understanding the mechanical behaviour of concrete's different phases and the developments in computational capability. Building on the recent advances, the work formulates 3D mesoscale models, reports and analyses the outcomes of the simulations. The paper is organized as follows: Section 2 describes the meso-structure geometric model and mesh generation; Section 3 presents the adopted models for material behaviour; Section 4 reports and discusses the main simulation results on the mechanical behaviour and on the transport behaviour before and after cracking; Section 5 summarises the conclusions from the work.

Meso-structure of concrete
An extensive investigation recently carried out by Wang et al. [1] shows that gravel aggregates have spherical/ellipsoidal shapes, while crushed aggregates have polyhedral shapes, and a number of voids can be easily observed in concrete at meso-scale from X-ray CT images (see [37]). Therefore, spherical/ellipsoidal aggregates are considered for gravel aggregates, while polyhedral aggregates are considered for crushed aggregates. The size distribution of aggregates described by the Fuller curve [38] is summarised in Table 1, and is used in this study [39]. It should be pointed out that only coarse aggregates larger than 2.36 mm are explicitly modelled. From Table 1, it can be observed that aggregates above 2.36 mm occupy 98.6%. The fine aggregates smaller than 2.36 mm are neglected in this study, and are considered together with cement matrix as mortar [41]. This keeps the model fidelity and avoids the computational difficulty caused by enormous number of elements. Voids with size ranging from 2 to 4 mm as presented in [42] are included in the concrete meso-structure.
The generation process is based on the Monte Carlo sampling method wherein the aggregates and voids are generated from prescribed distributions of their size, shape and volume fraction. A three-level hierarchical method is proposed to reduce the computational cost, and a computationally efficient and versatile in-house program heterogeneous material generator has been developed using MATLAB [43]. The detailed algorithms and procedure in generating both 2D and 3D numerical heterogeneous materials can be found in our previous work [1,44]. Concrete cubic samples with dimensions of 25 mm � 25 mm � 25 mm are generated. Fig. 1 shows the numerically generated models with gravel and crushed aggregates (aggregates in green, voids in red and ITZs in grey). Aggregate gradation in Table 1, aggregate volume fraction P agg ¼ 30%, void content P void ¼ 1%, aspect ratios for ellipsoidal aggregates and voids (the ratio of the major axis to the minor axis) , the minimum and maximum edges of the polyhedrons N min ¼ 8, N max ¼ 16, minimum space between the edge of an aggregate and the boundary of the concrete specimen γ 1 ¼ 0.5 mm and minimum gap width between any two aggregates γ 2 ¼ 0.5 mm, are adopted here [42].

Finite element mesh
The generated meso-structures of concrete, illustrated in Fig. 1, are meshed with the pre-processing functionality in commercial FE package ANSYS [45]. A purpose written in-house program was developed to identify the solid elements of ITZs between mortar and aggregates. The detailed meshing procedure is described as follows: Firstly, the volume representing aggregates and ITZs V a-I and aggregates only V agg , are created. The "overlapping" Boolean operation in the ANSYS pre-processor is applied to separate ITZs domain from aggregates and ITZs domain. Thus, the remaining volume is ITZs, denoted as V ITZ. Secondly, the volume representing voids V void and the entire concrete volume V all are created subsequently and the "overlapping" Boolean operation is utilized again to volume domains V ITZ , V agg, V void and V all, which result in separating V mor from V all .
The numerical models are meshed with solid hexahedral elements so that more realistic crack paths can be obtained. The ANSYS Parametric Design Language (APDL) programs in combination with ANSYS batch processing are developed to provide a powerful tool of automatic mesh generation. Finite element mesh (element size�0.5 mm) corresponding to the concrete with gravel/crushed aggregates and ITZs (aggregates in grey, mortar in blue and ITZs in purple) shown in Fig. 1 is generated and is shown in Fig. 2.

Mechanical behaviour of constituents
The coarse aggregates usually have a significantly higher strength than mortar. Therefore, a linear elastic material model is used to model the aggregates. For mortar and ITZs, concrete damage plasticity model in ABAQUS [46] is employed in the present study. The concrete damage plasticity model assumes that the concrete failure mechanisms are tensile cracking and compressive crushing which is characterised by damage plasticity (see Fig. 3). The cracking is controlled by two hardening variables, ε pl t and ε pl c , linked to failure mechanisms under tension and compression loading, respectively.
In this work, the analytical expression for normal concrete proposed by Guo and Zhang [47], and later adopted in Chinese Code GB50010 [48] is used for hardening and softening. For both mortar and ITZs under uniaxial tension, the stress-strain response follows a linear elastic relationship until the value of the failure stress σ t0 , which corresponds to the onset of micro-cracking in the concrete, is reached. Beyond the peak stress, the following equation is employed to describe the post-peak softening curve [48]: (1) where f t and ϵ t0 are the peak stress and peak strain, respectively, and α t is a coefficient calculated by α t ¼ 0.312f 2 t . For uniaxial compression, the stress-strain response is linear until the value of initial yield stress σ c0 is reached, and the plastic response is characterized by stress hardening followed by strain softening beyond the ultimate stress σ cu . The whole relationship is approximated by the following equations [48].
If σc where f c and ϵ cu are the peak stress and peak strain, respectively, and α a and α d are coefficients calculated by α a ¼ 2.4-0.0125f c and α d ¼ 0.157-f 0:785 c 0.905. It is assumed that the uniaxial stress-strain curves can be converted into stress versus plastic strain curves using the following equations, where ε pl t and ε pl c are equivalent plastic strain in tension and compres- Fig. 1. Numerically generated 3D models for concrete with gravel and crushed aggregates (examples with P agg ¼ 30% and P void ¼ 1%).
sion, ε ck t and ε in c are cracking strain and inelastic strain, and ε el 0t , ε el 0c , ε el t and ε el c are calculated using Equations (8) and (9), where d t and d c are damage variables in tension and compression, which evolves monotonically from 0, representing the undamaged material, to 1 which represents total loss of strength.
The stress-strain relations (σ-ε) under uniaxial tension and compression loading are given by Equations (10) and (11), respectively, where the strains in tension and compression can be represented as follows, In order to account for different effects under tension and compression loadings, two damage criteria are used: one for tension and a second for compression. Substituting Equations (5) and (6) into Equation (11) and using Equations (9) and (10) result in where b t and b c are constant factors extracted from experimental tests, Here the values 0.1 and 0.7 are selected for b t and b c , respectively, as suggested by Britel et al. [49].
The material properties used in this study are determined based on the data from relevant publications [6,12], and are listed in Table 2. In addition, five other parameters are needed for damage-plasticity model used in this study: Dilation angle: 35 � , flow potential eccentricity: 0.1, biaxial/uniaxial compression plastic strain ratio: 1.16, invariant stress ratio: 0.667, viscosity parameter: 0.0005 are adopted for both mortar and ITZs. Tensile and compressive stress data are provided as a tabular   function of strain, which is directly obtained from the constitutive relation of mortar and ITZs. The damage variables defined by Equations (12) and (13) are defined as a tabular function of cracking strain and inelastic strain for tension and compression, respectively.

Transport through constituents
In this work, concrete is considered as a multiple-porosity medium consisting of: (1) various-sized and disconnected voids, (2) cracks, (3) intermediate diffusivity ITZs, (4) impermeable and disconnected aggregates, (5) low diffusivity mortar matrix. The mortar matrix system contains a large number of small well-connected or isolated pores, whereas voids are larger pores which refer to the millimetre scale voids in the range of 2-4 mm in the mechanical analysis. The voids are connected to cracks directly or indirectly through mortar matrix. The interconnectivity between voids is through mortar matrix. The cracks pass through the voids and the interface between aggregates and mortar matrix.
The transport mechanisms of chloride in concrete are usually categorized into four groups [50]: diffusion driven by concentration difference, permeation driven by hydraulic pressure difference, migration driven by electrical potential difference and capillary suction due to moisture content (pressure) difference. Any one or combination of these mechanisms can govern the ingress of chloride ions into concrete. In general, diffusion is considered as the governing transport mechanism of chloride ions when concrete is fully saturated with water. In this paper, the concrete pore system is assumed to be fully saturated with water for the transport analysis after cracking. The porosities of mortar and ITZs, as well as voids volumes are assumed constant in time. The size and shape of fractures determined by the mechanical solution at a given load level do not change during the transport calculations, which means that the propagation and development of fractures are not affected by the fluid flow and only dependent on the mechanical effects. Based on these assumptions, Fick's second law is used to govern the diffusion process: where C is the chloride concentration [mol/m 3 ], D the chloride diffusion coefficient [m 2 /s], and x the spatial coordinate [m]. The diffusivity of uncracked cement matrix is 5.65 � 10 À 12 m 2 /s from the experiments [51], the diffusivity of aggregates is 0, and the diffusivity of ITZs is 5.65 � 10 À 11 m 2 /s which is in the range of measured values [51][52][53]. For diffusivity in cracks, experimental study [54][55][56] shows that chloride crack diffusion coefficient, D cr, depends on the crack width. Below a threshold crack width, cracks have little influence on transport, mainly because of the self-healing effects inside cracks. Above the threshold crack width, the chloride diffusion coefficient is equal to a constant value as the product of subsequent hydration cannot block the crack [57]. Different scholars have proposed different threshold values [54][55][56]. In this work, 80 μm is adopted for the crack threshold values according to the experimental data measured by Jang et al. [56].

Finite element analyses
3D concrete specimens with spherical aggregates and voids (L x � L y � L z ¼ 25 mm � 25 mm � 25 mm concrete cubic, P agg ¼ 30%, P void ¼ 1%, γ 1 ¼ γ 2 ¼ 0.5 mm) were first loaded under uniaxial compression and tension (see Fig. 4). The loading was displacement-controlled with zero prescribed vertical (X-direction) displacements at nodes on the bottom surfaces and non-zero vertical displacements at nodes on the top surfaces. Horizontal displacements for the same nodes are left free, except for the nodes at the right lower corner of the specimens, which are fixed to prevent rigid body translation.
The commercial software ABAQUS/Explicit was used for finite element analysis. The explicit method in ABAQUS was originally developed to solve dynamic problems, so it is important that the inertial forces do not affect mechanical responses and provide unrealistic dynamic results when performing quasi-static analysis. In principle, the loading time should be long enough to minimise any dynamic effects but not too long to cause excessive computational cost. After a trial and error analysis, loading displacement d ¼ 0.1 mm, loading time t ¼ 0.1s were found to be long enough to ensure the quasi-static loading condition [1,58]. The diffusion in sound and cracked concrete is analysed by ABAQUS as well. The interface elements are inserted along the surface of elements where the damage values d t /d c from mechanical analysis reach to 1 and having displacement discontinuity with neighbouring elements. These interface elements are used to represent the realistic cracking. The crack width is equal to the prescribed displacement from the mechanical analysis. In this work, new interface elements are proposed which can not only consider the diffusion in the normal direction but also in the tangential direction as shown in Fig. 5 compared to previous studies [59,60]. The interface element has zero-thickness with nine nodes.
Experimental study [56,[61][62][63] shows that chloride diffuses in crack with a different coefficient which is much larger than the coefficient in concrete matrix. The diffusion coefficient through a crack, D cr . varied with the crack width and had no relationship with the binder materials. The relation between D cr and crack width has been described by different artificially equations to fit the experimental data [56,[61][62][63]. The calculation of the crack diffusion coefficient, D cr [m 2 /s], is cited from Djerbi [61] in this work.

Macroscopic stress-strain behaviour
A typical 3D concrete specimen described in section 3 is modelled and loaded under uniaxial compression and tension. Fig. 6 (a) shows the predicted mean stress-strain curves in X direction for a 3D heterogeneous concrete subjected to uniaxial, monotonic compression/tension loading, typical experimentally observed response [64] and another numerical result [65]. The mean stress is calculated as the sum of reaction force of nodes on the bottom surface divided by the cross-section area of the specimen. Fig. 6 (b) shows a good correlation of stress-strain curve of a 3D heterogeneous concrete subjected to uniaxial, monotonic tension loading predicted in this study, with experimentally observed response [66] and the result from another numerical work [44]. They are qualitatively similar to each other with a clear peak and sharp initial post-peak drop followed by a long shallow tail. However, one should not compare the curves quantitatively as they describe specimens of different distribution, volume fraction of aggregates and voids, etc.
It can be seen that the stress-strain curve predicted in this study is in very good agreement with those from experimental test and the other numerical work. The curves show that typical concrete behaviour is characterised by an initial elastic region, a nonlinear hardening range before the peak stress, and a following softening branch. However, given the wide variation in the experimental data, and that the specimens are of different aggregate and void spatial distribution, calibration of the model to represent the results of a single experimental test or other numerical work is not productive. It should also be noted that the stressstrain curve predicted shows more brittle behaviour than the experimental result. This could be due to the delayed initiation of micro-cracks or increasing creep effects at high stress levels measured in tests. It is still unclear as there are no experiments in which micro-crack initiation and evolution in the pre-peak phase has been properly investigated.

Damage evolution and crack propagation
In order to illustrate the cracking process under compression, three special loading points, corresponding to one in the hardening range, one at peak stress and another at late softening range are investigated, the results are shown in left, middle and right columns in Fig. 7, respectively. The compressive and tensile damage of the specimen clearly shows the localization process leading to crack formation and growth during the loading. The different views of cracking patterns seen are attributed to the heterogeneous random distribution of aggregates and voids. The reason for the emergence of tensile damage, Fig. 7(b), in compression is the material heterogeneity with aggregates and voids located randomly. The vertical loading on the top and bottom of aggregates results in generating inclined compressive stress, which leads to the emergence of horizontal tensile stress between aggregates. It is worth noting that the diagonal compression failure predicted is associated with the compressive law where more brittle behaviour is assumed, while experimental tests sometimes show multiple vertical cracks before the development a single failure plane.
Due to the relative weakness of aggregate-mortar interfaces, compressive damage starts at the hardening stage, left column in Fig. 7 (a), and tensile damage starts at the peak stress in aggregate-mortar interfaces, middle column in Fig. 7(b). Around the peak stress, damage increases and micro-cracks (mostly the cracks in aggregate-mortar interfaces) start coalescing, middle column in Fig. 7. With the loading progressing into the softening region, both the compressive and tensile damage accelerate and the cracks propagate through mortar in the late softening stage, right column in Fig. 7. The predicted cracking evolution process is found to be in agreement with both experimental and numerical work in elucidating mechanisms of the micro-cracking and damage evolution in compression [65,67].
The detailed cracking patterns under uniaxial tension at peak stress, initial softening and late softening are presented in left, middle and right columns in Fig. 8, respectively. The tensile and compressive damage the specimen are used to clearly demonstrate the cracking initiation and propagation during the loading.
It can be seen that micro-cracks start emerging after reaching the peak stress, left column in Fig. 8, and tensile damage starts developing especially at the aggregate-mortar interfaces perpendicular to the loading direction, middle column in Fig. 8. As the loading displacement increases, some micro-cracks coalesce into dominant cracks, while the other cracks arrest due to stress redistribution, right column in Fig. 8. It should also be noted that a much smaller number of elements are under compressive damage in tensile loading, Fig. 8(b), compared to the elements under tensile damage for compression case, Fig. 7(b). This indicates that in compression both tensile and compressive behaviour of the concrete phases are activated, while in tension the concrete response is dominated by the tensile behaviour of different phases.

Effect of material heterogeneity
The effect of material heterogeneity is investigated by performing tension simulation on three concrete specimens with different aggregate   Fig. 9, the stress-strain curves are almost identical in the pre-peak range, however, the difference in peak strength is about 4% for tension, and larger differences are observed between the post-peak responses. Fig. 10 compares the damage and cracking patterns for the three different samples under tension. The different plots demonstrate that the distribution of aggregates and voids has a significant effect on the cracking patterns under uniaxial tension. Two typical macro-crack patterns for uniaxial tension as presented by image-based modelling [68] and cohesive zone modelling [44] are also observed in this study: Type I cracking with only one dominant crack, Fig. 10(a) and (b), and Type II cracking with two or more dominant cracks, Fig. 10(c). A larger load-carrying capacity in the softening range is observed for type II cracking, sample 3 in Fig. 10(b). This is explained by the emergence of multiple concurrently growing macro-cracks, leading to increased toughness.
Clearly, crack propagation and the emergent stress-strain response in concrete under tension are affected by its highly disordered mesostructure. The results obtained in this study indicate that conclusions based on a small number of analyses may be drawn for the pre-peak range; however, it should be taken carefully to draw conclusions for the peak strength and post-peak response, and the diffusivity of cracked concrete on the basis of a small number of analyses. Therefore, it is necessary to do a statistical analysis if the peak strength and softening process are what we are concerned. Further study will concentrate on developing parallel computing algorithms to reduce computational cost and improve efficiency for statistical analysis of damage and failure of a large number of 3D concrete specimens.

Effect of void content
The effect of void content on the fracture behaviour is investigated by fixing the positions of aggregates and aggregate volume fraction and varying the void content. All the other geometric and material parameters except the void content are fixed to the values used in Sections 4.1.1 and 4.1.2. As shown in Fig. 11, the effect of void content on the resulting stress-strain curves both under uniaxial compression and tension is pronounced, especially for peak and post-peak range. The loadcarrying capacity decreases as the void content increases for both compression and tension. This suggests that voids enable the cracks to propagate easily, acting as perforation unzipped by the cracks. The void distribution for three specimens with same aggregate distribution and aggregate volume fraction (P agg ¼ 30%), but different void content, i.e. 0%, 1% and 2% is shown in Figs. 12 and 13. It clearly shows that the voids, which as around the same size scale as aggregates exist in the concrete, should not be neglected when analysing the mechanical properties and fracture of concrete. It should also be noted that the distribution of voids may also have some effect on the resulting stressstrain curves, and a quantitative prediction could only be obtained by statistical analysis. Figs. 12 and 13 demonstrate further that the effect of void content on the cracking patterns is pronounced: the crack surface passes through larger number of voids as the void content increases. For both tension and compression, cracks tend to propagate along voids and neighbourhood of aggregates; preferable paths are chosen if given aggregate is relatively large. All observations suggest that it is important to consider voids, which have rarely been considered in meso-structural of concrete, when analysing the failure and damage behaviour.

Diffusion validation in sound concrete
In order to verify the transport model, the predicted effective chloride diffusion coefficients by the model in this work are compared with the measured values in Care [51]. The aggregate sizes and percentage where N is the total number of aggregates, r i is the radius of each aggregate, t is the ITZ thickness and V is the total volume of the model. In this case, for the same aggregate volume content but different aggregate size distribution, the larger aggregate surface area of aggregates the larger ITZ volume content will be obtained which has the same trend in Care [51].
The effects of aggregate volume fraction and ITZs volume fraction are studies as well in Fig. 14. The effective chloride diffusion coefficient decreases as the aggregate volume percentage increases for the same aggregate size distribution. While the effective chloride diffusion coefficient increases as the ITZs volume increases for the same aggregate volume fraction. The presence of aggregates attribute to the drop of the effective diffusion coefficient which induce a more tortuous path the transport of species, as well as a phase which is much more impermeable than the cement paste. While the formation of an ITZ around aggregates in mortar increases chloride ingress [69].
The difference between the simulated results and measured experimental values is due to that, firstly, the calculation of ITZ volume in Care [51] takes into account the overlap of ITZs while the calculation of ITZ volume in this work does not considered; and secondly, the chloride binding is neglected in the model.

Effect of ITZ diffusivity
The realistic thickness of ITZs in normal concrete may be different with the thickness used in this study. In order to study the effects of ITZs on the chloride ions diffusion and to avoid the computation difficulty caused by the element mesh of ITZs, the change of ITZs thickness can be reflected by the change of the diffusion coefficients of ions in ITZs. When  the aggregate size distribution is assigned as Section 2.1, the aggregates are considered as sphere and occupy 30% volume, and there are no voids in the model, the influence of ITZ diffusivity on chloride ions diffusion is analysed. The ITZ thickness was assumed to be constant and the diffusivity of ITZ is set as 0.1, 1, and 10 times of mortar matrix, respectively. Concentrations equal to 1 and 0 were applied to the inlet and outlet respectively. It should be noted here that the concentration gradient and assigned diffusivity values do not affect final results as they will be reported in terms of ratio of diffusivity of concrete to mortar matrix. This is due to the fact that the governing differential equations are linear and consequently the results are independent of these input values and adopted units. The simulated results are shown in Figs. 15 and 16.
The ITZ effect dominants the diffusion of chloride ions over the tortuosity effect with the increase of diffusivity of ITZ, as illustrated in Fig. 15. The contour lines of concentration in Fig. 16 become smoother and smoother when R (i.e. the ratio between diffusivity of ITZ and diffusivity of mortar) grows from 0.1 to 10. The iso-concentration lines in Fig. 16(b) nearly become flat lines when the diffusivity of ITZ is 10 times of mortar matrix. A considerable amount of chloride ions will diffuse through ITZ instead of mortar matrix. When the ratio, R, is equal to 1, the diffusivity of ITZ is the same as mortar matrix and the influence of ITZ on the chloride ions diffusion can be neglected. The concentration of chloride ions at the outer side of aggregate is higher than at in cement mortar at the same depth and the iso-concentration lines at the outer   side of two adjacent aggregate particles are outer convex lines. This is because that part of chloride ions will accumulate at the outer side of aggregate and other chloride ions will diffuse around aggregate through the mortar matrix alternatively. When the influence of ITZ is ignored, the iso-concentration lines are much more tortuous around aggregates than in case of R equal to 10. This is due to that the ITZ effect partly can counterbalance the tortuosity effect induced by aggregates.

Effect of void content
The effects of bigger voids on the concrete diffusivity before cracking are then analysed based on the same conditions in Section 4.1.4. Both the aggregates and voids morphology are considered to be sphere. The chloride ions diffusivities increase as the void content increases in Fig. 17.

Effects of crack width
From mechanical analysis in Section 4.1 and transport analysis in Section 4.2, the chloride ions diffusivities and cracking patterns vary as the change of parameters relevant to ITZ, aggregates and voids. Therefore, in this work, only the influence of crack width is studied in this work to verify the proposed model. The percentage and distribution of aggregates and voids are kept unchanged as Section 4.1.1. while only prescribed displacement changes to illustrate the effects on crack width on diffusion coefficients. The variation of the ratio D cracked /D sound with the crack width is illustrated in Fig. 18. D cracked is the diffusion coefficient in cracked specimen and D sound is the diffusion coefficient in sound specimen. D cracked /D sound increases as the crack width increases and is in the range of experimental values [61].

Conclusions
A multi-phase meso-scale model is proposed to study the sequential coupling of mechanical and transport analysis in concrete. From the results obtained the following conclusions can be drawn: (1) Mechanical part of the model is validated under both tension and compression by comparison of predicted and experimentally measured stress-strain curves (quantitatively) and of predicted and experimentally observed 3D cracking patterns (qualitatively). The results suggest that the model provides a good representation of the fundamental characteristics of concrete tensile and compressive responses. (2) Parametric analyses show that both aggregate and void content have pronounced effects on the macroscopic stress-strain response and failure patterns in compression and tension. Specifically, omitting voids from meso-scale models could lead to over prediction of tensile and/or compressive strength and unrealistic predictions for post-peak behavior and failure patterns. This in turn could lead to unrealistic basis for constructing cohesive laws for continuum modelling of concrete structures, should this be the purpose of meso-scale analysis.    and experimental measured effective chloride ions diffusion coefficients. The predicted effective diffusion coefficient of chlorides in sound concrete is very close to the measured experimental values under different aggregate volume fraction and size distribution. (4) The predicted effective diffusion coefficient of chlorides in sound concrete decreases as the aggregate volume increases for the same aggregate size distribution while it increases as the ITZs volume content and the diffusivity of ITZs increases. This is due to that the aggregates are more impermeable compared to other phases and they increase the tortuosity of transport path. The ITZs accelerates the transport as it has a larger diffusivity compared to mortar matrix. (5) The presence of larger voids has a greater effect on both mechanical response and transport properties. The load-carrying capacity decreases as the void content increases in both compression and tension. It is further shown that the void content has a large effect on the failure patterns. This leads to the conclusion that voids should be considered in meso-scale fracture modelling of concrete. Further study will concentrate on the behaviour of concrete under complex loading conditions, i.e. triaxial stress fields, present in different locations of engineering components, as means for deriving improved constitutive laws for continuum-based modelling of concrete structures. (6) The presence of microcracks has a greater effect on the diffusion.
The effective diffusion coefficient of chlorides in cracked concrete increases when the crack width above a threshold value increases. (7) The modeling approach described in this work is particularly useful for evaluating the effect of various parameters, such as the shape and size distribution of aggregate, volume fraction of ITZs, aggregates and voids, and microcracks, on mechanical response and transport properties, where laboratory experimentation alone would be difficult to isolate and quantify the effects. An important element of future research would be statistical analysis of the effects of shape, size, gradation, spatial distribution,   17. Effects of void content on chloride ions diffusivity. Fig. 18. Effects of crack width on chloride ions diffusivity (D cracked is the diffusion coefficient of cracked concrete while D sound is the diffusion coefficient of uncracked concrete; w is the crack width).
volume fraction of aggregate and void on the concrete behaviour, for which parallel computing algorithms are required.
Finally, it should be stressed that the present model does not take into account the effect of the overlap of ITZs and ionic binding. Thus, strictly speaking, improved model should be developed to predict the ionic diffusion coefficients more accurately.

Compliance with ethical standards
(in case of Funding) Funding: This study was funded by EPSRC via grant EP/N026136/1.

Declaration of competing interest
The authors declare that they have no conflict of interest.