P1Freely available online through the publisher-supported open access option.


Geological repository designs employ a multi-barrier approach. The materials, which include wasteforms, backfill and host rock, are typically porous quasi-brittle. Mechanical damage (e.g. nucleation and growth of microcracks) can result in significant changes in permeability. A knowledge of how the permeability is affected is critical to accurate modelling of radionuclide transport. This work proposes a novel 3D lattice-type model for the damage evolution in such materials, referred to as the site-bond model. Its advantages over previous models are that the shape of the lattice cell is physically realistic and that any macroscopic elastic response can be represented, including those of cementitious and geological materials. Damage accumulates as bonds fail upon reaching prescribed failure strengths. These are dictated by a predefined pore size distribution. Concrete is used as a study material. It is demonstrated that the model can predict the macroscopic stress–strain response under unconfined tension and compression with emergent non-linearity due to damage evolution. Ongoing work on the prediction of permeability changes with damage is discussed. This is based on the interaction between the model proposed here and a lattice model of the pore space.


The current plans for the long-term isolation of nuclear wastes involve the construction of a deep underground repository in suitable low-permeability rocks to ensure low mobility and high spatial dispersion of radioactive species at sufficient depths. In addition, several engineered barriers are considered to increase the retention of these species within the repository. These will be produced by encapsulating the wastes in cementitious or glassy forms, placing the waste forms in metallic containers, and subsequently backfilling the space between the containers and the host rock with cement or a clay-based material. The rock and all engineered barriers, except the containers, are porous brittle or quasi-brittle in nature. They contain one or several solid phases with different chemical compositions and mechanical properties along with a system of spatially distributed pores and/or other micro-defects. Their main function, to retain radionuclides, is to a large extent supported by the low permeability, which is dictated by the poor connectivity of the defect space.

Defect space changes, both in geometry and in connectivity, will result in permeability evolution. Geometry changes may result, for example, from electrochemical and/or bacterial processes eroding the accessible defect space. More detrimental, however, are the changes in connectivity. The main cause for these is the nucleation of additional defects, such as micro-cracks, and the growth and coalescence of defects. We will call this material damage. In brittle or quasi-brittle materials, unlike metals, the damage is initially distributed throughout the volume and its evolution depends on the distribution of these defects. Damage may evolve with changes in the mechanical loading and/or in the material mechanical properties. The former can be due to external processes, such as land motion or glaciation–melting cycles, or due to internal processes such as pressure build-ups from corrosion generated gases. The latter can be due to chemically or radiation-induced transformations of the solid phase. Although it is important to understand the damage-inflicting processes, it is equally important to be able to model the damage evolution of these materials and link this to the evolution of permeability. This is critical to the assessment of the long-time performance of the repository near-field.

Continuum modelling of distributed damage can be highly computationally demanding and is not possible without additional constitutive assumptions for the nucleation, growth, competition and coalescence of micro-defects. Particle-based methods, such as peridynamics (Silling, 2000) and smoothed particle applied mechanics (Hoover and Hoover, 2001), have been proposed to overcome these limitations. These, however, do not offer an evident way to distribute initial damage, e.g. in the form of pores, and link this to defect space models. This is because the methods are mainly concerned with the mechanical representation of crack nucleation and growth. More advantageous is the class of lattice-based models, developed specifically for quasi-brittle materials. Here, the solid is represented by a system of sites, connected appropriately by bonds that characterize relative deformations between solid-phase grains. Early works, based on 2D regular hexagonal lattices, provide a good illustration that the non-linearity in the macroscopic stress–strain curve for such materials is an emergent property from the underlying micro-cracking events (Schlangen and van Mier, 1992; Chang et al., 2002). However, the 2D image of damage is not credible, because the micro-crack growth and coalescence is a 3D phenomenon. A 3D model has been proposed based on a regular lattice with cubic cells (Schlangen, 2008). One problem with this lattice is that it is not physically realistic, neither in terms of the shape of the represented grains, or in terms of the shapes of the formed grain boundaries, nor in terms of the coordination with neighbouring grains. Moreover, this lattice cannot provide the macroscopic elastic response of the materials studied. In particular, this and other regular lattices with constant in length and shape bonds have been shown to represent isotropic elasticity only for materials with zero Poisson's ratio (Wang and Mora, 2008). Therefore, results obtained with such lattices could provide qualitative but mechanically unrealistic damage evolution. A lattice model, based on a cell known as the truncated octahedron, has been previously used for studies of intergranular stress corrosion cracking (Jivkov et al., 2006). This offers a more physically realistic description of solids at their meso-scale, the length scale dictated by the defect space, and can be tuned to yield any pre-defined macroscopic elastic behaviour, which is an essential prerequisite for the subsequent modelling of damage.

The aim of this work is to demonstrate the use of the novel lattice model for the prediction of damage evolution in porous quasi-brittle media. Concrete has been selected as the study material, because of its planned use in the repository, widely available data on the pore size distribution required for model construction, and data on its mechanical behaviour also required for model validation. In addition, we discuss the opportunities to link this evolution to the changes in permeability using a complementary lattice model of the defect space. The term site-bond model is introduced as a generic name for both the solid-phase damage model and the defect-space permeability model.

Model and method

The material meso-scale is considered to be a length scale relevant to a given mechanism of microscopic failure. The mechanism is related to the presence of micro-defects acting as failure initiators. A physical region of the material can be tessellated into cells, so that all micro-defects are distributed only at the interfaces between cells. For porous quasi-brittle materials the cells define solid-phase grains with pores or other defects residing at their boundaries. As we aim to construct a regular model of the microstructure, it is important to select the cells shape and coordination which are sufficiently representative of real materials. Monte Carlo studies with 3D Voronoi tessellations of space provide an insight into the selection of such a shape (Kumar, 1992). Of all regular solids that can fill the space compactly the best fit to the average shape found is the solid known as the truncated octahedron, bounded by six squares and eight regular hexagons. We use this as the unit cell of a regular lattice support (Fig. 1a) and construct a corresponding site-bond model (Fig. 1b) with sites at cell centres and bonds connecting sites of cells with common interfaces.

Note the presence of two type of bonds in the site-bond model, denoted by B1 and B2, with lengths L1 = L and 3L2, respectively, where L is the unit-cell size. The possibility to alter their stiffness coefficients independently makes the model sufficiently flexible to represent any predefined macroscopic elastic behaviour (Jivkov and Yates, 2012). Assuming local isotropy, the bonds are modelled with structural beam elements with circular cross sections with radii R1 and R2, respectively. The stiffness coefficients of the two types of bonds then depend on the ratios R1/L, R2/L, and Gb/Eb, where Eb and Gb are the moduli of elasticity and rigidity of the bonds.

We have selected plain concrete as a study material. Typical values for the macroscopic elastic and failure properties of this material are: compressive strength fc = 40 MPa; tensile strength ft = 3.3 MPa; modulus of elasticity E=5000fc; Poisson's ratio ν = 0.25 (Bortolotti, 2003). The linear elastic behaviour of this material is represented by a site-bond model with the following bond properties: R1/L = 0.1; R2/L = 0.35; Gb/Eb = 0.61; and Eb = 66292 MPa. Pores are distributed hypothetically on the boundaries between cells. This corresponds to defects in bonds, which can cause bond failure in separation (fracture) or sliding (shear) mode. Proximity to failure is measured by a failure factor, F. Separation can occur for bonds in tension, for which F = σ /σc, where σ > 0 is the total stress acting on the pore and σc is the critical stress that depends on the pore size. Sliding can occur for bonds in compression, for which F = (|τ| + σ tanϕ) /τc, where τ and σ <0 are the shear and the compressive stresses acting on the pore, ϕ is the angle of internal friction, and τc is the critical shear stress that depends on the pore size. The typical value for concrete ϕ = 37° is used here. A bond fails if F >1. We have used published data for the size distribution of pores (Winslow and Liu, 1990) and have converted it to strength distribution assuming Griffiths failure criterion, i.e. σc and ~1r, where r is the pore size. This lead to the following Weibull-type strength distributions: σc = Sc [−ln(1 – p)]1/m and τc = Tc [−ln(1–p)]1/m, with shape m = 2.73. The scale parameters used are Sc = 6.25 MPa and Tc/Sc = 4; the latter represents the typical ratio between mode II and mode I fracture toughness for concrete. The bond strengths are distributed using a generator for uniformly distributed random numbers 0 ⩽ p < 1.

We use a model occupying the region (20L, 20L, 20L) with respect to a coordinate system (X1, X2, X3) oriented along the principal axes of the unit cell. The deformation/damage behaviour of the model is studied as a sequence of bond-failure events. At each simulation step i, the model is subjected to a fixed tensile or compressive macroscopic strain along one of the axis, εm, representing unconfined tension or compression, respectively. The corresponding macroscopic stress, σmi, is determined from the reaction forces. The bond with maximum failure factor Fi >1 at this step is removed from the system for the next step. The macroscopic stress and strain at bond failure are σi=σmiFi and ei=εmFi, respectively, because of the linear relation between bonds and model behaviour. This is a linear version of the return-mapping algorithms used in plasticity modelling. The sequence (εi, σi) defines the macroscopic stress–strain behaviour of the material. The macroscopic modulus of elasticity evolves as Ei = σii, and a damage parameter can be defined in the standard way as Di = 1 – Ei/E, where E = σ00 = 31623 MPa is the initial modulus.


The macroscopic stress–strain behaviour of the concrete under tensile (assumed positive) and compressive (assumed negative) loading predicted by the model is shown in Fig. 2. The value of the macroscopic damage parameter, D, is depicted at several locations along the curves. The damage under tension and compression in terms of micro-crack population at D = 0.05 is illustrated in Fig. 3. The micro-cracks are shown as discs. Those formed by separation are given in yellow and those formed by sliding are given in blue. The entire populations are embedded in cubes with sides parallel to the model.

The dominant failure mechanism under tension is separation, where bonds parallel (B1) or inclined (B2) to the load fail in tension. The macroscopic damage evolves quickly with the population count. Under compression, a mixture of separation failures, where bonds perpendicular (B1) or inclined (B2) to the load fail in tension, and sliding failures, where bonds inclined (B2) to the load fail in shear, is observed. The macroscopic damage evolves slowly with the population count. A comparison of the two evolutions is shown in Fig. 4a, where the damage for the two cases is plotted against the fraction of failed bonds. For the compressive load, the sliding and separation failures causing damage are shown in Fig. 4b as fractions of the total number of failures.


The results presented in Fig. 2 are very encouraging. We have used simply the elastic modulus and the Poisson's ratio of plain concrete to calibrate the elastic properties of the bonds and an experimental pore-size distribution for this material to obtain distributions of failure strengths of these bonds. These seem to be sufficient to make a very good prediction of the tensile and the compressive stress–strain behaviour of the material. In tension: (1) the experimentally observed non-linearity prior to failure is clearly reproducible; (2) the predicted tensile strength (peak of curve) is very close to the experimental value; and (3) the experimentally observed post-failure softening behaviour is also captured. This behaviour is almost exclusively due to local separations normal or inclined to the applied load. This results in a rapid increase of macroscopic damage with the number of local failures, as seen in Figs 3a and 4a. In compression: (1) the almost linear behaviour prior to failure is reproduced; (2) the non-linearity increases within a progressive damage region; and (3) the predicted compressive strength (peak of curve) is slightly larger than the experimental one, but if the stress at the onset of progressive damage is taken as compressive strength than this is very close to the experimental strength. This behaviour is due to a mixture of dominating local separations parallel or inclined to the applied load and an increasing number of sliding failures in compressed inclined bonds (Fig. 4b). This results in a slow increase of macroscopic damage with the number of local failures, as seen in Figs 3b and 4a.

The problem formulated and solved for illustration is based on minimal microstructure data and the simplest possible macroscopic representation of the material as a homogeneous isotropic continuum. The model, however, allows for a number of additional features to be readily introduced including: (1) anisotropy by assigning different elastic or failure properties to bonds in different directions (seven distinct directions can be used). (2) heterogeneity by assigning different elastic or failure properties to the bonds within different domains of the model; this can be used to represent material grains of different sizes. (3) material texture by assigning different cell sizes in the principal directions of the truncated octahedron. (4) initial damage by removal of a set of bonds from the assembly, either randomly or according to a particular rule. Currently, the model is being developed and tested on concrete with in-house obtained pore space data, from X-ray CT, and under various complex loadings, including extension, compression and shear meridians following initial confining pressures of up to the concrete unconfined compressive strength. The results will be reported in a future publication.

The main advantage of the proposed model is that it is conceptually simple and easy to implement, yet sufficiently powerful to include a wide range of phenomena in a wide class of materials. Furthermore it offers an elegant way to determine the effect of damage evolution on the permeability. This can be accomplished with a secondary defect-space model that supplements and interacts with the primary solid-phase model. In the secondary model the sites represent micro-defects, located in the middle of the primary model bonds. The bonds represent links between micro-defects, such as throats between pores or nucleated micro-cracks. There is no requirement that all sites and bonds in the defect-space model are present initially. This will depend on the size distribution of micro-defects for the particular problem. The secondary model is a site-bond type model, conceptually similar to the pore network models used in petrology (Blunt et al., 2002). The permeability can be determined as a solution of the defect-space model (Jivkov and Olele, 2011). The link between the primary and the secondary model is a subject of an ongoing work. This will allow for all failure events in the solid-phase model to be accounted for correspondingly in the defect-space model so that the evolution of permeability with damage can be derived.


We have proposed a regular lattice model that conceptualizes porous quasi-brittle materials at their meso-scale in which: (1) the unit cell represents the average shape that can be found when tessellating such materials into defect-free domains; and (2) defects reside between cells and their size distribution determines the distribution of local strengths within the model components. We have demonstrated that the model can be successfully used for studies of damage evolution due to mechanical loading, presently unconfined tension and compression. This evolution involves the natural creation of micro-cracks, their growth, interaction and coalescence. From a micro-structural perspective, the model can be used for defining anisotropic and heterogeneous domains, both in terms of the solid-phase and the defect-space, thus being able to represent a wide class of porous or fractured brittle or quasi-brittle media. From a material performance perspective, the model can be used for obtaining realistic damage evolution, hence integrity assessment, as well as permeability evolution with the help of a secondary defect-space model.