Hydro-mechanical network modelling of particulate composites
Ignatios Athanasiadis Simon J. Wheeler Peter Grassl
cond-mat.mtrl-sci



    





  Ignatios Athanasiadis, Simon J. Wheeler and Peter Grassl^*

  School of Engineering, University of Glasgow, Glasgow, UK

  ^*Corresponding author: Email: peter.grassl@glasgow.ac.uk, Phone: +44 141 330 5208

 
  



Keywords: Microcracking, particulate composites, aggregate restrained shrinkage, mass transport, lattice, hydro-mechanical, network model, periodic boundary conditions



§ ABSTRACT

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 <cit.> and has shown to increase mass transport properties such as permeability and sorptivity <cit.>.
Numerically, the initiation of microcracking due to particle restrained shrinkage was studied in <cit.>.
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 <cit.>.
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 <cit.> 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 specimen 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 <cit.>, 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 <cit.>.
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 <cit.>.



§ 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 <cit.>.
The transport constitutive model used was Darcy's law combined with a cubic law <cit.> 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 <cit.> 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 freedom 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 <cit.>.
However, the method in <cit.> 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 <ref> 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 𝐱 and 𝐱' are the coordinate vector of the original point and one of the image points, respectively. Furthermore, 𝐌 is the translation matrix defined as 𝐌 = diag [1 + k_ x a, 1 + k_ y b, 1 + k_ z c] where k_ x, k_ y, k_ z∈{-1,0,1}.
The k_ x, k_ y and k_ z coefficients define the direction of the shift from the original point to the image point.
The coordinates of the 26 image points are the result of the coordinate translations in (<ref>) for all k_ x, k_ y, k_ z combinations except for the case where k_ x = k_ y = k_ z = 0.
In Figure <ref>, the cell with two of its 26 neighbours is shown.
The points I and J are examples of two randomly placed points satisfying the minimum distance requirements.
Points I' and J' are one of 26 sets of periodic image points of I and J, respectively, whereby I' was generated by a translation 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 <cit.>.
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 <cit.>.
In Figure <ref>a, a Delaunay tetrahedron and the Voronoi facet associated with Delaunay edge i-j are shown.
The structural elements were placed on the Delaunay edges with their mid-cross-sections defined by the facets of the Voronoi polyhedra (Figure <ref>b).
Analogous to the structural network, the transport elements were placed on the edges of the Voronoi polyhedra, with their cross-sections formed by the facets of the Delaunay tetrahedra (Figure <ref>(c)).



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 <ref>, an example of two elements crossing the cell boundaries is shown.
Here, the intersections of the elements I-J' and I'-J with the boundaries of the periodic cell are presented with a circle and a cross in the plane of the boundary.
The parts of the elements that lie outside the cell boundaries are shown with dashed lines and those inside by solid lines.
In the present periodic cell approach, all degrees of freedom of nodes inside the cell boundaries were involved in the system of equations that was used for determining the unknown degrees of freedom.
For elements crossing the boundary, the degrees of freedom of the nodes outside the cell were determined from those of the periodic image of the node inside the cell and additional information in the form of average strain and pressure gradients for the structural and transport problem, respectively. 

Examples of structural and transport networks generated from the same set of random points are presented in Figures <ref> and <ref>, respectively.





For the structural network, the structural elements are shown in Figure <ref>a and the edges of the mid-cross-sections of the structural elements in Figure <ref>b.

The transport elements of the transport network are presented in Figure <ref>a and the edges of the transport mid-cross-sections are shown in Figure <ref>b.
From these examples of structural and transport networks, several interesting features of this coupled periodic cell are visible.
Firstly, both structural and transport elements (Figures <ref>a and <ref>a) are either located inside the cell or cross its boundaries.
No elements are entirely placed outside the periodic cell. 
The edges of the cross-sections of structural and transport elements are either inside the cell, cross the boundary or lie entirely outside the periodic cell.
Therefore, the network of structural elements in Figure <ref>a is smaller than the network of edges of the cross-sections of transport elements shown in Figure <ref>b. Analogue to the structural network, the network of transport elements in Figure <ref>a is smaller than the network of the edges of cross-sections of the structural network in Figure <ref>b.
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 <cit.>.
Details regarding how the conductivity was calculated are presented in section <ref>. 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 <cit.>, which is 


    ∇σ^ c = 0


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 <cit.>.
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 (<ref>) for the structural element shown in Figure <ref>(b) is


    𝐊𝐮_ e = 𝐟_ s


where 𝐊 is the stiffness matrix, 𝐮_ e are the vector of degrees of freedom and 𝐟_ 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 𝐮_ e = {𝐮_ t^T, 𝐮_ r^T }^T, where 𝐮_ t = {𝐮_ ti^T, 𝐮_ tj^T}^T = { u_ xi, u_ yi, u_ zi, u_ xj, u_ yj, u_ zj}^T and 𝐮_ r = {𝐮_ ri^T, 𝐮_ rj^T}^T = {ϕ_ xi, ϕ_ yi, ϕ_ zi, ϕ_ xj, ϕ_ yj, ϕ_ zj}^T. 
These degrees of freedom 𝐮_ t and 𝐮_ r are used to determine displacement discontinuities 𝐮_ C = {u_ Cn, u_ Cp, u_ Cq}^T at point C by rigid body kinematics <cit.> as 


    𝐮_ C = 𝐁𝐮_ e = 𝐁_ 1𝐮_ t + 𝐁_ 2𝐮_ r


where 𝐁 = {𝐁_1, 𝐁_2}^T, 𝐁_1 and 𝐁_2 are two matrices containing the rigid body information for the nodal translations and rotations, respectively, which are


    𝐁_ 1 =
    [ -𝐈  𝐈 ]


and


    𝐁_ 2 =
    [     0 -e_ q  e_ p     0  e_ q -e_ p;  e_ q     0  -h/2 -e_ q     0  -h/2; -e_ p   h/2     0  e_ p   h/2     0 ]


Here, 𝐈 is a 3×3 unity matrix.

In (<ref>), 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 <ref>b).
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. 

The displacement jump 𝐮_ C in (<ref>) is transformed into strains ε = {ε_ n, ε_ p, ε_ q}^T = 𝐮_ C/h, where h is the length of the structural element.
The strains are related to stresses σ = {σ_ n, σ_ p, σ_ q}^T by means of a material stiffness 𝐃 = (1-ω) 𝐃_ e, where 𝐃_ e = diag{E, E, E }. Here, E is the Young's modulus and ω is the damage variable, which is further discussed in Section <ref>.
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


    𝐊 = Ah[ 𝐁_ 1^ T𝐃𝐁_ 1 𝐁_ 1^ T𝐃𝐁_ 2; 𝐁_ 2^ T𝐃𝐁_ 1 𝐁_ 2^ T𝐃𝐁_ 2 ]
      +
      [               0               0;               0 𝐁_ 1^ T𝐊_ 𝐫𝐁_ 1 ]


Here, 𝐊_ r is a matrix containing the rotational stiffness at point C defined as


    𝐊_ r = (1-ω)Eh[ I_ p    0    0;    0  I_1    0;    0    0  I_2 ]


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 (<ref>) 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 <cit.>.

The above element formulation for structural elements entirely located in the periodic cell is identical to the one described in <cit.>.
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 𝐄 = {E_ x, E_ y, E_ z, E_ yz, E_ zx, E_ yx}^T. 
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 <ref>.
The translations of a node outside the cell is


    u'_ x = u_ x + ak_ xE_ x + ck_ zE_zx + bk_ yE_yx



    u'_ y = u_ y + bk_ yE_ y + ck_ zE_yz



    u'_ z = u_ z + ck_ zE_ z

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 <ref>. Node J' is outside the cell and its periodic image J is inside the cell.
Making use of (<ref>), (<ref>) and (<ref>) 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


    [  𝐮_I; 𝐮_J';  𝐫_I; 𝐫_J' ]
    =
    𝐓_ m[ 𝐮_I; 𝐮_J; 𝐫_I; 𝐫_J;   𝐄 ]


where 𝐮_I, 𝐫_I, 𝐮_J, 𝐫_J, 𝐮_J' and 𝐫_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 𝐓_ m is of size 12 × 18 and has the form

    𝐓_ m
    =
    [     𝐈     0     0     0     0     0;     0     𝐈     0     0 𝐤_ 21 𝐤_ 22;     0     0     𝐈     0     0     0;     0     0     0     𝐈     0     0;       ]

The sub-matrices 𝐤_ 21 and 𝐤_ 22 are 3 × 3 matrices which contain information about the transformation of the nodal translations due to the average strains.
They are defined as

    𝐤_21
     =
    [ ak_ x     0     0;     0 bk_ y     0;     0     0 ck_ z;       ]

and

    𝐤_22
     =
    [     0 ck_ z bk_ y; ck_ z     0     0;     0     0     0;       ]


If (<ref>) is combined with (<ref>) for calculating the displacement jump, the transformation matrix 𝐓_ m multiplies matrix 𝐁 from the right.
It follows from duality that the internal forces must be multiplied by 𝐓_  m^ T from the left, before the evaluation of the equilibrium conditions.
Hence, the original 12 × 12 stiffness matrix 𝐊 of the non-periodic element is now transformed into the 18 × 18 matrix 𝐓_ m^ T𝐊𝐓_ m.

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 <cit.>, 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 σ = {σ_ n, σ_ p, σ_ q}^T as 

    σ = (1-ω) 𝐃_ e(ε-ε_ p - ε_ s) = (1-ω) σ̅

where ω is the damage variable, 𝐃_ 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 (<ref>), flow rule (<ref>), evolution law for the hardening variable (<ref>) and loading unloading conditions (<ref>):


    f = F(σ̅,κ)



    ε̇_ p = λ̇∂ g∂σ̅



    κ̇ = λ̇ h_κ



    f ≤ 0 λ̇≥ 0 λ̇ f = 0


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 = √(σ̅_ s^2 + σ̅_ t^2) is 

    f = {[        α^2σ̅_ n^2 + 2 α^2(f_ c - αβ f_ t)1+αβ q σ̅_ n + σ̅_ q^2 - 2 α^2 f_ cf_ t + α^2 (1-αβ) f_ t^21+αβ q^2                                                                                                             ; σ̅_ n^2β^2 + 2 f_ c - αβ f_ t1+αβ q σ̅_ n + σ̅_ q^2 + (1-α^2 β^2) f_ c^2 -2 αβ(1+αβ) f_ c f_ tβ^2 (1+αβ) q^2                                                                                                              ].

where f_ t and f_ c are the tensile and compressive strengths, respectively, and α and β are the friction angles shown in Fig. <ref> for f = 0 and q=1, which controls the hardening.



It is defined as

    q = exp(κA_ h)

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 (<ref>) are the same as those of the yield surface f except that α is replaced by ψ:


    g = {[ ψ^2 σ̅_ n^2 + 2  ψ^2 (f_ c - ψβ f_ t)1+ψβ q σ̅_ n + σ̅_ q^2                                                            ;         σ̅_ n^2β^2 + 2 f_ c - ψβ f_ t1+ψβ q σ̅_ n + σ̅_ q^2                                                             ].

The smaller ψ is, the smaller is the ratio of normal and shear components of plastic strains for σ̅_ n≥f_ c - ψβ f_ t1+ψβ q.
The function h_κ in the evolution law in (<ref>) is chosen as 

    h_κ =  |∂ g∂σ_ n |

which is the absolute value of the normal component of the direction of the plastic flow.

The damage variable in (<ref>) is determined by means of the damage history variable

    κ_ d = ⟨ε_ pn⟩

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

    𝐰_ c = h ( ε_ p + ω(ε_ n - ε_ pn))

For pure tension, the crack opening simplifies to

    w_ c = h ( ε_ pn + ω(ε_ n - ε_ pn))
 
where h is the length of the network element (Figure <ref>).
The stress-crack opening curve is

    σ_ n = f_ texp(-w_ cw_ f)

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 (<ref>) equal to the first component of (<ref>), a nonlinear equation of the damage ω is obtained, which is solved using the Newton-Raphson method.
For modelling the dependence of transport properties on cracking, permeability, which is part of the transport model described in section <ref>, is made dependent on the absolute value of the crack opening  

    w̃_ c = |w_ c|


The structural constitutive model requires eight input parameters.
The Young's modulus of the lattice material E controls the macroscopic Young's modulus.
The parameters of the plasticity part are f_ t, f_ c, α, β, ψ and A_ h.
Finally, G_ F controls the amount of energy dissipated during cracking.



 §.§ Transport network

For the transport part of the model, a 3D network of 1D transport elements is used to discretise the stationary transport equation <cit.>


    div(k grad P_ f) = 0


Here, P_ f is the fluid pressure, k is the conductivity. In (<ref>), 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 (<ref>) for a 1D transport element shown in Figure <ref>c is


    𝐤_ eP_f = 𝐟_ e


where 𝐤_ e is the one-dimensional element conductivity and 𝐟_ e are the nodal flow rate vector <cit.>. The degrees of freedom of the transport elements are the fluid pressures P_ f = {P_ f1, P_ f2}^T.
Within the context of a one-dimensional finite element formulation <cit.>, Galerkin's method is used to construct the elemental conductivity matrix as 

    𝐤_ e = k A_ th_ t[  1 -1; -1  1 ]

Here, A_ t is the  mid-cross-sectional area and h_ t the length of the transport element shown in Figure <ref>c.



For elements, which cross the faces of the cell, the element formulation is explained based on Figure <ref> used earlier for the structural element. For instance, for the element I'J in Figure <ref>, 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

    [ P_f I';   P_fJ ]
    =
    𝐓_ t[     P_fI;     P_fJ; Δ P_fx/a; Δ P_fy/b; Δ P_fz/c ]

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 <ref>.
Furthermore, 𝐓_ t is a transformation matrix of size 2 × 5 and has the form

    𝐓_ t =
    [     1     0 ak_ x bk_ y ck_ z;     0     1     0     0     0;       ]

For combining (<ref>) with (<ref>), the transformation matrix 𝐓_ t multiplies matrix 𝐤_e from the right.
It follows from duality that the internal flux must be multiplied by 𝐓_ t^ T from the left, before the evaluation of the balance condition.
The conductivity matrix is evaluated as 𝐓_ t^ T𝐤_ e𝐓_ t, where the original conductivity matrix 𝐤_ 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 


    k = k_0 + k_ c


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

    k_0 = ρκ_0μ

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 (<ref>) models the increase of conductivity due to cracking using a cubic law based on the concept of flow through parallel plates <cit.> with a reduction factor ξ for the presence of roughness of the wall surface <cit.>.
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 <cit.>. However, since this is an important part of the model of the present study, it is shown here once more.
The term k_ c is


    k_ c = ξρ12 μ A_ t∑_i=1^3w̃_ci^3 l_ci


where w̃_ci and l_ci are the equivalent crack openings and crack lengths (see Figure <ref>) 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 𝐰_ c defined in (<ref>). 
The relation in (<ref>) expresses the well known cubic law, which has shown to produce good results for transport in fractured geomaterials <cit.>. 
The way how crack openings in the structural elements influence the conductivity of a transport element is schematically shown in Figure <ref>.
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 (<ref>) in proportion to their equivalent crack widths and the crack lengths. This crack length (shown by blue double lines in Figure <ref>) is defined as the length from the midpoint of the structural element to the centroid C_ t of the transport element cross-section.



§ ANALYSES



 §.§ Introduction

The coupled structural transport network model described in section <ref> 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 <ref>.
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 <ref> according to <cit.>.





Three groups of analyses were carried out.
Firstly, the structural response of the matrix material subjected to mechanical loading was analysed by means of direct tension and a hydrostatic compression test using the structural periodic cell.
With these analyses key factors of the performance of the structural constitutive model were demonstrated.
This study was required since the model used here differs from the one in the previous two-dimensional study in <cit.> 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 <ref>.

For direct tension, the periodic cell was subjected to monotonically increasing average axial strain E_ x introduced in section <ref>. 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 <ref>a and b, respectively. 



The structural network approach exhibits the typical response of cohesive-frictional quasibrittle materials subjected to tension in the form of softening, i.e. decreasing stress with increasing displacement (Figure <ref>a), and localised deformation (Figure <ref>b).
The peak stress is larger than the tensile strength input f_ t, because for an irregular network arrangement the individual elements are subjected to a combination of normal and shear stresses. For the present input, the shear strength is greater than the tensile strength. For shear at zero normal stress at the onset of hardening (called here f_ q), the elastic limit is f_ q = 2f_ t for the values of the parameters α and β in Table <ref>.

For hydrostatic compression, the normalised average stress σ_ v/f_ c versus strain response and the crack patterns at the stage marked in the stress-strain curve are shown in Figure <ref>a and b, respectively. Here, σ_ v = (σ_ x + σ_ y + σ_ z)/3 and ε_ v = (ε_ x + ε_ y + ε_ z)/3. 


The overall response in hydrostatic compression is initially elastic followed by elasto-plastic hardening.
The deviation from the elastic response occurs at σ_ v = -f_ c, which corresponds to the onset of the yielding in the constitutive model for negative normal stress only (σ_ n<0 and σ_ q = 0).
The cracks patterns, using the definition of the crack opening in (<ref>), are distributed within the periodic cell without showing any patterns of localised deformations, which is typical for quasi-brittle materials subjected to hydrostatic compression.
The crack opening consists only of plastic strains, since no damage occurs in these hydrostatic compression analyses. 

The two examples of direct tension and hydrostatic compression demonstrate that the structural constitutive model is capable of describing the response in tension and compression realistically.
This is important for the particle restrained shrinkage analyses in sections <ref> and <ref>, in which complex tensile and compressive stress states play an important role.
The responses for both tension and compression are insensitive to the element size.
For direct tension, the irregularity of the network affects the results, because the inelastic displacements are localised in an element size dependent region.
However, global result in the form of the stress-displacement curve is insensitive to the element size.
For hydrostatic compression, the results converge with network refinement.
For large element sizes, the stiffness is overestimated because elements crossing the boundary contribute stronger to the stiffness of the cell.
Since in the analyses in sections <ref> the network size is varied, this insensitivity to the network size is important. 



 §.§ 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 (a=b=c in Figure <ref>) 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 

    a = (n_ pπ6ρ_ p)^1/3 d

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 <ref>.
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

    κ_ yy = Q_ yμA_ yρΔ P_ y/L_ y

where A_ y = a× b is the cross-sectional area in y-direction (Figure <ref>).
The permeability component in (<ref>), which is one of nine components of the matrix of permeability <cit.>, was used here to assess the influence of particle size on permeability for microcracking induced by particle restrained shrinkage.

In Figure <ref>, the mean of the permeability κ_ yy normalised by the intrinsic permeability of the matrix κ_0^m 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 <ref>.
Once cracking has been initiated (ε_ s<-0.1 %), permeability is strongly influenced by the 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 <ref>a 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 (<ref>) 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 <cit.>, 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 <cit.>.
This particle size dependent crack opening results in a strong dependence of conductivity on particle size, since the cubic law in (<ref>) 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 <ref>b 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 <ref>b, it can be seen that the flow paths are orientated preferentially in the vertical direction.
Two of the crack planes in Figure <ref>a 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 <ref> is in agreement with the two-dimensional results in <cit.>.



 §.§ 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 <ref>).
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 <ref>.
As in the single particle analyses in section <ref>, 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 (<ref>).
The mean of the permeability of ten random analyses versus the shrinkage strain is presented in Figure <ref> 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 <ref>.
However, the permeability at the final stage is similar for the random and regular particle arrangements.
In Figure <ref>, the crack patterns of one of the thick specimens (c=75 mm) is shown for three stages of applied shrinkage strain (ε_ s = -0.25, -0.375 and -0.5 %).
The blue spheres indicate the position of the particles.
The yellow polygons show mid-cross-sections of elements in which the crack opening increases at this stage of analysis and is greater than 10 μm.
Colours refer to the online version.
The spacing between the crack planes is determined by the size of the particles and their spatial arrangement.
For all three stages of applied shrinkage strains, the crack network connects randomly placed particles.
Very similar observations were made in the two-dimensional study in <cit.>.
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 <ref>.
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 <ref>a, none of the transport elements exhibits flow greater than Q_0.
The crack openings of the crack networks in Figure <ref>a 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 <ref>b, a network of transport elements exceeding the threshold is visible.
These transport elements are located on the crack planes shown in Figure <ref>b.
Not all transport elements on crack planes exhibit high flow, because of the variation of the crack openings of individual elements.
In the final stage, the majority of transport elements located on crack planes shown in Figure <ref>b conduct flow greater than Q_0. 
At this stage, the crack network is fully developed and the cracks have opened up so much that transport elements along the crack planes provide the majority of the flow through the specimen.


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 <ref> 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 <ref>, 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 <ref>.
For the network with the smallest specimen thickness (c=25 mm) in Figure <ref>a, the number of elements with high flow is smaller than for the other two thicknesses in Figure <ref>b and c. This appears to be in contradiction to the results in Figure <ref>, 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 <ref> and <ref> 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 <ref> 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 <cit.>. This influence of nonuniform shrinkage strain on microcracking and conductivity changes will be investigated in future studies.




§ ACKNOWLEDGEMENTS

The numerical analyses were performed with the nonlinear analyses program OOFEM <cit.> extended by the present authors.
The authors acknowledge funding received from the UK Engineering and Physical Sciences Research Council (EPSRC) under grant EP/I036427/1 and funding from Radioactive Waste Management Limited (RWM) (http://www.nda.gov.uk/rwm), a wholly-owned subsidiary of the Nuclear Decommissioning Authority. RWM is committed to the open publication of such work in peer reviewed literature, and welcomes e-feedback to rwmdfeedback@nda.gov.uk.

elsarticle-harv


equationsection

