Hydro-mechanical network modelling of particulate composites

Differential shrinkage in particulate quasi-brittle materials causes microcracking which reduces durability in these materials by increasing their mass transport properties. A hydro-mechanical three-dimensional periodic network approach was used to investigate the influence of particle and specimen size on the specimen permeability. The particulate quasi-brittle materials studied here consist of stiff elastic particles, and a softer matrix and interfacial transition zones between matrix and particles exhibiting nonlinear material responses. An incrementally applied uniform eigenstrain, along with a damage-plasticity constitutive model, are used to describe the shrinkage and cracking processes of the matrix and interfacial transition zones. The results showed that increasing particle diameter at constant volume fraction increases the crack widths and, therefore, permeability, which confirms previously obtained 2D modelling results. Furthermore, it was demonstrated that specimen thickness has, in comparison to the influence of particle size, a small influence on permeability increase due to microcracking.


Introduction
Microcracking due to particle restrained shrinkage significantly increases the permeability of porous quasi-brittle materials, which often reduces the durability of these materials.
For instance for cementitious composites, microcracking due to particle (aggregate) restrained shrinkage has been experimentally observed in [3,25,16,26,17] and has shown to increase mass transport properties such as permeability and sorptivity [25,26]. Numerically, the initiation of microcracking due to particle restrained shrinkage was studied in [8,13,10]. In some of these studies, it was shown that the width of cracks produced by particle restrained shrinkage depends strongly on the size of particles [8,10]. In the two-dimensional numerical modelling, particles were often idealised as cylindrical particles and all cracks were assumed to penetrate completely the specimen thickness. Therefore, in these two-dimensional analyses, increase of crack width resulted in a significant increase in permeability.
In experiments of irregular particulate composites such as concrete, permeability was measured by applying a unidirectional pressure gradient of either water or gas across the specimen thickness. In these thicker specimens, it can be assumed that complicated 3D fracture networks are generated due to differential shrinkage. It appears to be reasonable to expect that part of these networks are not connecting the opposite sides of the specimen and, therefore, do not equally contribute to the increase of permeability. Even for crack paths connecting opposite sides of the specimen, the crack widths along the path might vary. Therefore, it can be expected that the thicker the specimen is, the smaller the increase of permeability due to cracking will be. For instance, recent experimental studies in [26] for nonuniform drying shrinkage showed that permeability depends on specimen thickness. This dependence on specimen thickness could be due to the nonuniformity of the shrinkage strain, resulting in thickness dependent patterns of microcracking in the specimen or due to the variation of crack openings along random crack planes.
The aim of the present study was to investigate numerically the separate influences of spec-imen thickness and particle size on the increase of permeability due to particle restrained shrinkage induced microcracking. For this purpose, a new coupled hydro-mechanical periodic network approach consisting of coupled structural and transport networks was developed. Periodic cells are known in the area of homogenisation [19,11], where it has been shown that periodic boundaries result in faster convergence of properties with increasing cell size than boundaries subjected to displacement or traction conditions. One of the new features of the present periodic approach is that not only the periodic displacement/pressure-gradient conditions were applied, but also the three-dimensional network structure has been chosen to be periodic for both the structural and the transport networks. Peviously, combinations of periodic network structure and periodic displacement conditions were developed only for two-dimensional structural networks [7]. The new three-dimensional coupled periodic approach allows for describing fracture patterns and resulting increases in conductivity independent of boundaries. For the transport network, the constitutive models were based on Darcy's law combined with a cubic law to model the influence of cracking on permeability [5].

Method
The present numerical method for investigating the influence of particle restrained shrinkage induced microcracking on transport properties was based on three-dimensional periodic hydro-mechanical networks of structural and transport elements. For the constitutive model of structural elements, a damage-plasticity model was used [6]. The transport constitutive model used was Darcy's law combined with a cubic law [24] to model the increase of permeability due to fracture. The new feature of the present method is the extension of the coupled network approach proposed in [5] to a periodic cell of the shape of a rectangular cuboid, which uses periodic network structures and periodicity requirements for the displacements, rotations and pressures. In addition to nodal degrees of freedom in the form of displacements, rotations and fluid pressures of nodes inside the cell, average strain and pressure gradient components are used to solve for unknown degrees of free-dom of nodes outside the cell. The use of a periodic cell with periodic network structure and periodicity requirements for the degrees of freedom has the advantage that for the mechanical network, crack patterns are independent of the cell boundaries, which is not the case for boundary value problems with either displacements or tractions applied to the boundaries. The formulation of the hydro-mechanical periodic cell approach is conceptually based on work reported in [7]. However, the method in [7] was limited to a two-dimensional structural network. Here, a three-dimensional coupled hydro-mechanical periodic network approach was proposed.

Discretisation
The periodic dual network approach was based on Delaunay and Voronoi tessellations of a set of points placed randomly within a rectangular cuboid shown in Figure 1 with thick lines. The points were placed sequentially while enforcing a minimum distance d min between all placed points. Trial points that fail the minimum distance criterion were rejected. The placement was terminated once the number of trials for placing one point exceeds the limit N iter .
Once the placement of points was completed, 26 periodic image points were generated for each successfully placed point within the cell using the translation rule where x and x are the coordinate vector of the original point and one of the image points, respectively. Furthermore, M is the translation matrix defined as with k x = −1 and k y = k z = 0, and J with k x = 1 and k y = k z = 0.
All points within the cell and all periodic image points are used for the Delaunay and Voronoi tessellations. The Delaunay tessellation decomposes the domain into tetrahedra whose vertices coincide with the randomly placed points. The Voronoi tessellation divides the domain into polyhedra associated with the random points [21]. Each polyhedron is the subset of the domain in which points are closer to the placed point that is associated with the polyhedron than all the other placed points. Facets of Voronoi polyhedra form subsets of the 3D space, in which every location is equidistant from a pair of placed points and nearer to these two points than to any other point. The edges of Delaunay tetrahedra connect pairs of placed points of Voronoi polyhedra with common facets.
Delaunay and Voronoi tessellations were used to define the structural and transport elements [5]. In Figure 2a Because of the periodic image points, the tessellated space is larger than the main cell.
Therefore, edges of Delaunay tetrahedra and Voronoi polyhedra cross the cell boundaries.
In  Figure 3b. For the coupling of the two networks, cross-section edges which are entirely outside the cell require special consideration. In the present network approach, a one way coupling approach was used in which crack openings obtained from the structural network were used to compute the conductivities of transport elements. Therefore, for the transport network, crack openings associated with cross-sectional edges outside the cell were assumed to be equal to those of the corresponding cross-sectional edges located inside the cell [2]. Details regarding how the conductivity was calculated are presented in section 2.3.2. The only input parameters required for the discretisation of the periodic cell are the minimum distance d min and the maximum number of trials to place one point N iter . These parameters control the average lengths of structural and transport elements.
The greater N iter is, the smaller is the standard deviation of the element lengths up to the stage at which the domain is almost saturated with points and and increase of N iter will result in small changes of the number of placed points.

Structural network
The three-dimensional structural network was designed to approximate the quasi-static equilibrium equation without body force [23], which is where ∇ is the divergence operator and σ c is the continuum stress.

Structural element
The structural element formulation for elements which lie entirely in the periodic cell is identical to the one presented in [5]. However, for elements crossing the boundary of the periodic cell, a new element formulation was introduced. To be able to explain this new feature, the standard formulation is presented first. The discrete version of (2) for the structural element shown in Figure 2(b) is where K is the stiffness matrix, u e are the vector of degrees of freedom and f s are the acting forces. The formulation of the structural element is presented in the local coordinate system, i.e. the coordinate system (x, y and z) of the nodal degrees of freedom coincides with the coordinate system (n, p and q) of the quantities used for evaluating the constitutive response. Each node has 3 translational (u x , u y and u z ) and 3 rotational (φ x , φ y and φ z ) degrees of freedom. The degrees of freedom of a structural element with nodes i and j are grouped in translational and rotational parts as These degrees of freedom u t and u r are used to determine displacement discontinuities u C = {u Cn , u Cp , u Cq } T at point C by rigid body kinematics [12] as where B = {B 1 , B 2 } T , B 1 and B 2 are two matrices containing the rigid body information for the nodal translations and rotations, respectively, which are and Here, I is a 3 × 3 unity matrix. In (6), e p and e q are the eccentricities between the midpoint of the network element and the centroid C in the directions p and q of the local coordinate system, respectively ( Figure 2b). The local coordinate system is defined by the direction n, which is parallel to the axis of the element, and p and q, which are chosen as the two principal axes of the mid-cross-section.
where h is the length of the structural element. The strains are related to stresses σ = Here, E is the Young's modulus and ω is the damage variable, which is further discussed in Section 2.2.2. For the present elastic stiffness matrix D e , Poisson's ratio equal to zero is obtained and the structural network is elastically homogeneous under uniform modes of straining.
For the case that the global coordinate system coincides with the local one, the element stiffness matrix is Here, K r is a matrix containing the rotational stiffness at point C defined as Here, I p is the polar moment of area, and I 1 and I 2 are the two principal second moments of area of the cross-section. The factor (1 − ω) in (8) ensures that the rotational stiffness reduces to zero for a fully damaged cross-section (ω = 1). The stiffness matrix is then expressed in the global coordinate system by means of rotation matrices as described for instance in [18].
The above element formulation for structural elements entirely located in the periodic cell is identical to the one described in [5]. For elements crossing the cell boundaries, a special formulation is required. For these elements, the degrees of freedom of the nodes outside the cell are determined from the degrees of freedom of the periodic image inside the cell and the average strain Here, E x , E y and E z are the average normal strains in the x, y and z direction, respectively and E yz , E zx , E yx are the average engineering shear strain components. For an illustration of the coordinate system x, y and z, see Figure 1. The translations of a node outside the cell is where the translation presented without and with the prime symbol are those of the nodes located within and outside the cell, respectively. Note that the contributions of the average shear strains E zx and E yx have been included only in the displacements in the x direction and the contribution of E yz has been included only in the displacement in the y-direction. This is justified because rigid body rotations of the entire cell are arbitrary. One node of the network is fully fixed in order to prevent rigid body rotation and translation.
Consider the element IJ in Figure 1. Node J is outside the cell and its periodic image J is inside the cell. Making use of (9), (10) and (11) and assuming that φ xJ = φ xJ , φ yJ = φ yJ and φ zJ = φ zJ , the transformation rule giving the translations and rotations of the two ends I and J of a structural element IJ crossing a cell boundary is where u I , r I , u J , r J , u J and r J are the vectors containing translational and rotational degrees of freedom of nodes I, J and J , respectively. The node J is the periodic image of point J inside the cell. The transformation matrix T m is of size 12 × 18 and has the form The sub-matrices k 21 and k 22 are 3 × 3 matrices which contain information about the transformation of the nodal translations due to the average strains. They are defined as and If (12) is combined with (4) for calculating the displacement jump, the transformation matrix T m multiplies matrix B from the right. It follows from duality that the internal forces must be multiplied by T T m from the left, before the evaluation of the equilibrium conditions. Hence, the original 12 × 12 stiffness matrix K of the non-periodic element is With the present approach, average stresses and strains are prescribed by means of the additional six average strain components introduced in the formulation of the periodic cell. Finally, the total number of degrees of freedom are six times the number of nodes positioned within the cell boundaries plus six additional degrees of freedom, which correspond to the six average strain components.

Structural material
The constitutive model for the structural material is based on a damage-plasticity framework [6], which is capable of reproducing the important features of the response of quasibrittle materials in tension and compression. The strains are related to the nominal stress where ω is the damage variable, D e is the elastic stiffness, ε p = ε p n , ε p p , ε p q T is the plastic strain andσ is the effective stress. Furthermore, ε s = {ε s , 0, 0} T is the shrinkage strain which was used in this study to initiate microcracking.
The plasticity model used to determine the effective stress is independent of damage.
The model is described by the yield function (17), flow rule (18), evolution law for the hardening variable (19) and loading unloading conditions (20): Here, f is the yield function, κ is the hardening variable, g is the plastic potential, h κ is the evolution law for the hardening parameter andλ is the rate of the plastic multiplier.
The yield function of the two stress variablesσ n andσ q = σ 2 s +σ 2 t is where f t and f c are the tensile and compressive strengths, respectively, and α and β are the friction angles shown in Fig. 5 for f = 0 and q = 1, which controls the hardening. It is defined as where A h is an input parameter. For the onset of plastic flow κ = 0 and q = 1.
The stress dependent parts of the plastic potential g in the non-associated flow rule in (18) are the same as those of the yield surface f except that α is replaced by ψ: The smaller ψ is, the smaller is the ratio of normal and shear components of plastic strains The function h κ in the evolution law in (19) is chosen as which is the absolute value of the normal component of the direction of the plastic flow.
The damage variable in (16) is determined by means of the damage history variable where . denotes the McAuley brackets (positive part of operator). The function of the damage variable is derived from the stress-crack opening curve in pure tension (σ n > 0, σ q = 0). For the damage-plasticity constitutive model, the vector of crack opening components is defined as For pure tension, the crack opening simplifies to where h is the length of the network element ( Figure 1). The stress-crack opening curve is where w f controls the initial slope of the exponential softening curve. It is related to the area under the stress-crack opening curve G F as w f = G F /f t . Setting (28)

Transport network
For the transport part of the model, a 3D network of 1D transport elements is used to discretise the stationary transport equation [15] div (kgradP f ) = 0 (30) Here, P f is the fluid pressure, k is the conductivity. In (30), gravitational effects are neglected.

Transport elements
Analogue to the structural network, different element formulations are used for elements which are located entirely in the cell and those that cross one of the faces of the cell.
For those inside the cell, the discrete form of (30) for a 1D transport element shown in Figure 6: Influence of cracking on transport. Figure 2c is where k e is the one-dimensional element conductivity and f e are the nodal flow rate vector [14,4]. The degrees of freedom of the transport elements are the fluid pressures Within the context of a one-dimensional finite element formulation [14], Galerkin's method is used to construct the elemental conductivity matrix as Here, A t is the mid-cross-sectional area and h t the length of the transport element shown in Figure 2c.
For elements, which cross the faces of the cell, the element formulation is explained based on Figure 1 used earlier for the structural element. For instance, for the element I J in Figure 1, the node I is outside the cell and J is inside the cell. The periodic image of I is I which is located inside the cell, which is used together with average fluid pressure gradients to determine the fluid pressure of the node outside the cell. The nodal fluid pressures of an element crossing the cell boundary are where ∆P fx /a, ∆P fy /b and ∆P fz /c are the average fluid pressure gradients along the x, y and z directions respectively. Here, a, b and c are the dimensions of the periodic cell shown in Figure 1. Furthermore, T t is a transformation matrix of size 2 × 5 and has the form For combining (34) with (31), the transformation matrix T t multiplies matrix k e from the right. It follows from duality that the internal flux must be multiplied by T T t from the left, before the evaluation of the balance condition. The conductivity matrix is evaluated as T T t k e T t , where the original conductivity matrix k e is transformed from a 2 × 2 matrix to a 5 × 5 one. The global conductivity matrix is assembled normally except for six rows and columns that relate the global degrees of freedom to its conjugate reaction flow rates.
As a result, the periodic cell can be subjected to arbitrary combinations of average flux or gradients. The total number of unknown degrees of freedom is the total number of the nodes located in the interior of the periodic cell plus three global degrees of freedom controlling the average flux or pressure gradient in the three directions.

Transport materials
The conductivity matrix for the material of the transport elements is where k 0 is the initial conductivity of the undamaged material and k c is the change of conductivity due to fracture. The conducticity of the undamaged material is where ρ is the density and µ is the dynamic viscosity of the fluid, and κ 0 is the intrinsic permeability. In this work, the density and dynamic viscosity was set to ρ = 1000 kg/m 3 and µ = 0.001 Pa s, respectively, which corresponds to the values commonly used for water.
The second term k c in (35) models the increase of conductivity due to cracking using a cubic law based on the concept of flow through parallel plates [24] with a reduction factor ξ for the presence of roughness of the wall surface [1]. A detailed description of the definition of k c and its dependence on the crack openings of the structural network, which was used in this study, has already been presented in [5]. However, since this is an important part of the model of the present study, it is shown here once more. The term wherew ci and l ci are the equivalent crack openings and crack lengths (see Figure 6) of neighbouring structural elements, which are located on the edges of the cross-section, and ξ is a reduction factor which considers the reduction of flow for cracks with rough surfaces compared to that between smooth parallel plates.
Here,w c is the magnitude of the crack opening w c defined in (27). The relation in (37) expresses the well known cubic law, which has shown to produce good results for transport in fractured geomaterials [24]. The way how crack openings in the structural elements influence the conductivity of a transport element is schematically shown in Figure 6. For instance, for the transport element o − p, three structural elements (i − k, k − j and i − j) bound the cross-section of the transport element. Thus, the conductivity will be influenced by these three elements according to (37) in proportion to their equivalent crack widths and the crack lengths. This crack length (shown by blue double lines in Figure 6) is defined as the length from the midpoint of the structural element to the centroid C t of the transport element cross-section.

Introduction
The coupled structural transport network model described in section 2 was applied to the analysis of particle restrained shrinkage of prisms made of a particulate quasi-brittle material consisting of particles, matrix and interfacial transition zones. The matrix and interfacial transition zones are made of permeable quasi-brittle cohesive-frictional material, where the interfacial transition zone is weaker and more permeable than the matrix.
The particles are elastic, and stiffer and less permeable than the matrix. The input parameters for the three material phases are presented in Table 1. The properties of the different phases are mapped onto a coupled periodic network. Network elements with both nodes within the same particle are given the property of particles. Those elements with nodes in different particles or one node in a particle and the other in the matrix are assigned the properties of the interfacial transition zone. Here, it is assumed that the thickness of the interfacial transition zone is much smaller than the length of the element, so that the Young's modulus of the elements crossing the interfacial transition zone is determined as the harmonic mean of those of matrix and particle. The strength of the element is determined by the strength of the interfacial transition zone. Finally, for elements with both nodes outside particles, the matrix material properties are used.
For all elements, the same structural and transport constitutive models were used with input parameters shown in Table 1 according to [8].
Three groups of analyses were carried out. Firstly, the structural response of the matrix  This study was required since the model used here differs from the one in the previous two-dimensional study in [8] and it is important that the numerical method provides network-independent results in tension and compression. In the second group, periodic cells with a centrally located single particle were analysed. Particle restrained shrinkage was modelled by subjecting matrix and interfacial transition zones to uniform incrementally increasing eigenstrain, while keeping the force resultants of the entire specimen at zero. The particle, which was not subjected to eigenstrain, restrained the matrix and interfacial transition zones. Therefore, this process is called particle restrained shrinkage.
After every increment of applied eigenstrain, the permeability of the specimens was evaluated by a stationary transport analysis with a fluid pressure gradient across the specimen applied in one direction. In these single particle analyses, the particle diameter was varied at constant particle volume fraction to investigate the influence of particle diameter on changes of permeability due to microcracking induced by restrained shrinkage. The last group consists of coupled analyses of specimens with multiple randomly placed particles of constant diameter and volume fraction for varying specimen thickness in the direction in which the pressure gradient is applied. With these analyses, the influence of specimen thickness on changes of permeability due to microcracking induced by particle restrained shrinkage were investigated. In the following sections, the three groups of analyses are discussed in detail.

Uniaxial tension and hydrostatic compression
Before investigating microcracking induced by particle restrained shrinkage, the performance of the structural network approach was investigated by means of direct tension and hydrostatic compression analyses. For this, a cubic cell with edge length of a = b = c = 5 cm was discretised with the periodic structural network approach using three network sizes of d min = 8, 4 and 2 mm. The number of iterations for the network generation was N iter = 10000. The material parameters for all network elements in this cube were the one of the matrix phase shown in Table 1.
For direct tension, the periodic cell was subjected to monotonically increasing average axial strain E x introduced in section 2.2.1. The average stress components corresponding to the other average strain components were set to be zero. The resulting normalised stress-displacement curve and crack pattern for the fine network at the stage marked in the stress-displacement curve are shown in Figure 7a and b, respectively. The structural network approach exhibits the typical response of cohesive-frictional quasibrittle materi-  Table 1.

Size of particles
In the second part of the study, the influence of size of inclusions on increase of permeability due to cracking was analysed. For these analyses, a coupled periodic cubic cell Figure 1) with a single centrally arranged particle was used. The volume fraction was kept constant at ρ p = 0.137, while varying the particle size as d = 4, 8, 16 mm. For the volume of a spherical particle V p = πd 3 /6 and the cell volume V cell = a 3 , the volume fraction for n p particles is ρ p = (n p V p ) /V cell . Combining these expressions and solving for the specimen length gives Therefore, for constant volume fraction ρ p , the size of the periodic cell decreases with decreasing particle size. For instance, for n p = 1 and d = 16 mm, the specimen length results in a = 25 mm. For the discretisation of the network, the ratio of the size of the cell and minimum distance was chosen as a/d min = 12.5 for all particle sizes investigated.
Thus, the average network element length decreases with decreasing particle size. This change of element size should not affect the results strongly, as it was shown in section 3.2.
For all analyses N iter = 10000 was used.
For the coupled analyses, a uniform shrinkage strain of ε s = −0.5% was applied in 100 increments to matrix and interfacial transition zone. The particles were not subjected to the shrinkage strain. For each particle size, ten analyses with random network generations were performed. After every increment of shrinkage strain, the permeability was determined by applying a unit fluid pressure gradient in the y-direction (∆P y /L y = 1).
Here, ∆P y and L y = a = b = c are the fluid pressure difference and the length in the y-direction, respectively. The total flow Q y in the y-direction resulting from this fluid pressure gradient was used to determine the macroscopic permeability component of the cell in the y-direction as where A y = a × b is the cross-sectional area in y-direction ( Figure 1). The permeability component in (39), which is one of nine components of the matrix of permeability [20], was used here to assess the influence of particle size on permeability for microcracking induced by particle restrained shrinkage.
In Figure 9, the mean of the permeability κ yy normalised by the intrinsic permeability of the matrix κ m 0 versus shrinkage strain ε s of ten random analyses is shown. The areas next to the curves show plus/minus one standard deviation. The scatter originates from the irregularity of the structural and transport networks. At early stages of the analyses (ε s > −0.1 %), microcracking due to particle restraint does not occur, so that there is no visible influence of particles size on permeability on the log-scale used in Figure 9. Once cracking has been initiated (ε s < −0.1 %), permeability is strongly influenced by the  Figure 9: Particle size: Permeability κ yy normalised by the intrinsic permeability of the matrix κ m 0 versus shrinkage strain ε s for three particle diameters at constant particle volume fraction. particle size at constant particle volume fraction. The greater the particle is, the greater is the increase of permeability. This strong dependence of permeability on particle size at constant volume fraction is explained by the crack patterns which are generated by the shrinkage of matrix and interfacial transition zone and the restraint that the particle provides. Crack patterns for the largest particle size (d = 16 mm) are shown in Figure 10a at the final increment of shrinkage strain. In this figure, cracks are visualised by yellow polygons representing mid-cross-sections of elements in which the equivalent crack width w c defined in (29) is greater than 10 µm. Colours refer to the online version. For all analyses of periodic cells with a single particle size, the shrinkage strain applied to matrix and interfacial transition zone results in overall regular localised crack patterns with three distinct crack planes aligned with the directions of the Cartesian coordinate system used for the periodicity of the cell. In an earlier two-dimensional study in [8], these type of regular crack patterns were obtained by placing multiple layers of inclusions in a regular pattern in a bigger specimen. Here, because a cell with periodic boundary conditions is used, these characteristic crack patterns for a regular inclusion arrangement are obtained for a cell with a single inclusion. The crack openings depend strongly on particle size. The greater the particle size is, the greater is the volume of material associated with this particle, which will be subjected to the eigenstrain and which will crack if the eigenstrain is sufficiently large. Therefore, the greater the particle size, the greater is the crack opening and the smaller is the crack length. This reasoning is in agreement with the two-dimensional numerical results reported in [8]. This particle size dependent crack opening results in a strong dependence of conductivity on particle size, since the cubic law in (37) was used to related crack opening to permeability due to cracking. Therefore, the increase in crack openings dominates the decrease in crack length.
In addition to the crack patterns, the main flow paths through the crack network is visualised in Figure 10b by showing transport elements in blue in which the flow is greater than the threshold Q 0 = 2.5 × 10 −16 kg/s. This value is equal to the entire flow through a cell of the same cross-section made of undamaged matrix material and subjected to a unit pressure gradient. In Figure 10b, it can be seen that the flow paths are orientated preferentially in the vertical direction. Two of the crack planes in Figure 10a are orientated so that they result in an increase of the mass transport in the y-direction. Most of the transport elements exceeding the threshold are located on these two crack planes. The strong dependence of permeability on particle size in Figure 9 is in agreement with the two-dimensional results in [8].

Specimen thickness
In the third part, the influence of specimen thickness on random particle arrangements in the direction of the pressure gradient was investigated for microcracking induced by particle restrained shrinkage. The aim of this part of the study was to investigate if potentially disconnected crack networks result in reduction of permeability with increasing specimen thickness. The cross-section of the rectangular periodic cell was chosen as a = b = 50 cm and the specimen thickness was varied as c = 25, 50 and 75 mm ( Figure 1). Furthermore, particle diameter and density were chosen as d = 16 mm and ρ p = 0.137, respectively. The minimum distance for the background network was chosen as d min = 2 mm. For each specimen thickness, ten analyses with random periodic particle and network arrangements were carried out using the coupled network approach presented in section 2. As in the single particle analyses in section 3.3, a shrinkage strain of ε s = −0.5 % was applied incrementally to matrix and interfacial transition zone elements. After every increment, the permeability component κ yy was determined as explained in (39). The mean of the permeability of ten random analyses versus the shrinkage strain is presented in Figure 11 in the form of lines with symbols for the three specimen thicknesses. Note that the mean results for c = 50 and 75 mm are almost indistinguishable. The coloured areas around the mean curves represent plus/minus one standard deviation. The permeability increases strongly with increasing shrinkage strain, as it was already observed for the single particle analyses. The strong increase of permeability is the result of the creation of random crack networks. These are more complex than the regular patterns obtained in section 3.3. However, the permeability at the final stage is similar for the random and regular particle arrangements. In Figure 12, the crack patterns of one of the thick specimens (c = 75 mm) is shown for three stages of applied shrinkage strain (ε s = −0. 25 Figure 11: Specimen thickness: Permeability κ yy normalised by the intrinsic permeability of the matrix κ m 0 versus shrinkage strain ε s for three specimen thicknesses at constant particle volume fraction and particle diameter.  [8]. The greater the shrinkage strain, the denser is the crack network and the greater are the crack openings. In addition, the flow network is shown in Figure 13. Blue lines indicate transport elements in which the flow is greater than Q 0 , which is the threshold used earlier for the single particle analyses. At the first stage in Figure 13a, none of the transport elements exhibits flow greater than Q 0 . The crack openings of the crack networks in Figure 12a are, while already localised, not large enough to increase the flow through the transport element sufficiently to exceed the threshold Q 0 . For the second stage in Figure 12b, a network of transport elements exceeding the threshold is visible. These transport elements are located on the crack planes shown in Figure 12b. Not all transport elements on crack The specimen thickness has overall only a small influence on the increase of permeability compared to the influence of the particle size and shrinkage strain. The mean permeability for c = 50 and 75 mm are almost identical. Only for the specimen with c = 25 mm, a greater permeability than for c = 50 and 75 mm was obtained. This difference is smaller than the standard deviation of the analyses with c = 25 mm. The greater the thickness is, the smaller is the standard deviation of the permeability, because some of the irregularities of the crack planes are averaged out along the specimen thickness. Crack patterns for specimens of three different thicknesses are shown in Figure 14 for the final stage (ε s = −0.5). For all three thicknesses, crack patterns connecting the particles are fully formed. Qualitatively, there is little difference between these crack patterns. In Figure 15, the network of transport elements in which the flow is greater than the threshold Q 0 is shown for a shrinkage strain of ε s = −0.5. For all three specimens, the network of transport elements exceeding this threshold have a very similar structure. The high flow elements are positioned on the crack planes as shown in Figure 14. For the network with the smallest specimen thickness (c = 25 mm) in Figure 15a, the number of elements with high flow is smaller than for the other two thicknesses in Figure 15b and c.
This appears to be in contradiction to the results in Figure 11, in which the specimen with the smallest thickness exhibits the greatest increase in conductivity due to cracking.
Nevertheless, the crack and flow patterns in Figures 14 and 15 are only one of ten random particle arrangements. These illustrations are useful for assessing the type of crack and flow patterns, but cannot be used for a quantitative comparison analyses with different thicknesses. For such a comparison, the mean values shown in Figure 11 should be used. In addition, for the smallest specimen thickness, the highest standard deviation was reported.

Conclusions
A three-dimensional coupled structural transport network model was used for the analysis of microcracking in heterogeneous materials due to particle restrained shrinkage.
A new coupled periodic network approach has been proposed, which combines periodic displacement/pressure-gradient conditions with periodic structural and transport networks. This new approach was applied to investigate the influence of particle size and specimen thickness on permeability for particle restrained shrinkage. The results of the three-dimensional analyses show that a change of particle size at constant volume fraction has a very strong influence on the increase of permeability due to microcracking for the case of particle restrained shrinkage. The influence of a change specimen thickness on increase of permeability for a pressure gradient in the direction of the changed specimen thickness was shown to be small in comparison.
This study was limited to shrinkage of matrix and ITZ applied uniformly across the periodic cells. In many applications, the shrinkage strain will be nonuniform through the specimen because of for instance a humidty gradient away from the specimen surface [9].
This influence of nonuniform shrinkage strain on microcracking and conductivity changes will be investigated in future studies.