Engineering Fracture Mechanics xxx (2013) xxx-xxx

Contents lists available at SciVerse ScienceDirect Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech M

ELSEVIER

Pore space and brittle damage evolution in concrete

Andrey P. Jivkova,c'*, Dirk L. Engelberg b,c, Robert Stein d,e, Mihail Petkovskie

a School of MACE, The University of Manchester, Manchester M13 9PL, UK b School of Materials, The University of Manchester, Manchester M13 9PL, UK

c Research Centre for Radwaste and Decommissioning, The University of Manchester, Manchester M13 9PL, UK d Nuclear FiRST Doctoral Training Centre, Universities of Sheffield and Manchester, Sheffield SI 3JD, UK e Civil and Structural Engineering, The University of Sheffield, Sheffield SI 3JD, UK

ARTICLE INFO ABSTRACT

A novel lattice model is proposed for linking experimentally measured porosity of concrete to damage evolution and the emergent macroscopic behaviour. Pore sizes are resolved by X-ray CT and distributed at lattice bonds. The mechanical behaviour of bonds is elastic-brittle with failure criterion dependent on local forces and pore sizes. Bond failures provide the only non-linear effect on the macroscopic response. Results are compared to several experimental load cases. They show good agreement of stress-strain response at lower stress levels and expected differences at peak stresses. The framework allows for future development of models with plasticity and time-dependent effects.

© 2013 Elsevier Ltd. All rights reserved.

Article history: Available online xxxx

Keywords: Concrete Lattice model Porosity

X-ray tomography Brittle fracture

1. Introduction

This work is part of an ongoing research programme on the performance of cement-based materials for nuclear plant and radwaste application. In some cases such materials will have predominantly radionuclide retaining function (e.g. waste-forms, backfills) and in others predominantly structural function (e.g. Advanced Gas-cooled Reactor pressure vessel, excavation support structures). For example, a critical role of a repository for radioactive waste is to ensure minimal release of radioactive species to the geosphere over very long times. An engineering safety case would demand fundamental understanding of the long-term evolution of the macroscopic properties of these materials. Key to the retaining function of the repository is the evolution of the transport properties, such as permeability and diffusivity. These are dictated by the 3D pore space in the materials used, which is characterised by the sizes and the connectivity of the pores present (see e.g. [1]). Consequently, the evolution of the macroscopic transport properties is governed by changes in the pore space.

The ongoing programme aims at developing, and validating experimentally, predictive models for the evolution of transport properties with pore space changes. This entails a modelling approach based on a practical, sufficiently realistic and modifiable 3D pore space representation. These requirements can be met to a large extent by the so called pore network models, where the pore space is described by a system of pores with various sizes some of which are connected by throats with various sizes [2,3]. Changes in the pore space may result from chemical, electrochemical or bacterial effects as well as from mechanical damage, such as microcracking. When the pore space changing mechanism is defined, the evolution of the macroscopic transport properties can be evaluated by linking the pore network model of the pore space to appropriate model for the selected mechanism [4]. This work focuses on the mechanically induced microcracking as a pore space changing mechanism in concrete. We seek a microstructure-informed 3D model, where the pore space is explicitly represented in

* Corresponding author at: Research Centre for Radwaste and Decommissioning, The University of Manchester, Manchester M13 9PL, UK. Tel.: +44 1613063765; fax: +44 1612003723.

E-mail address: andrey.jivkov@manchester.ac.uk (A.P. Jivkov).

0013-7944/$ - see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.engfracmech.2013.05.007

A.P.Jivkov et al./Engineering Fracture Mechanics xxx (2013) xxx-xxx

Nomenclature

Li, L2

Ri R2 S

v Vb rf sf

rmax smax £i ri ro

pore radius

material length scale (characteristic distance for pore failure) confining pressure

random number uniformly distributed macroscopic modulus of elasticity modulus of elasticity of lattice beams nodal forces, i = 1, 2, 3

cumulative probability function for pore sizes

macroscopic bulk modulus of elasticity

lattice spacing in principal directions

lengths of principal and octahedral beams in lattice

bending moment in lattice beam

critical (failure) bending moment

normal force in lattice beam

critical (failure) normal force

radius of beam with circular cross sections

radii of principal and octahedral beams with circular cross sections

shear force in lattice beam

critical (failure) shear force

twisting moment in lattice beam

critical (failure) twisting moment

coordinate axes, i = 1, 2, 3

nodal displacements, i = 1, 2, 3

shear to normal strength ratio

macroscopic Poisson's ratio

Poisson's ratio of lattice beams

critical (failure) normal stress in bond

critical (failure) shear stress in bond

maximum normal stress in beam due to bending

maximum shear stress in beam due to torsion

principal strains, i = 1, 2, 3

principal stresses, i = 1, 2, 3

"ideal" tensile strength (failure stress of material without defects) non-dimensional pore size effect parameter

terms of experimentally determined pore size distribution, which can be used to inform a pore network model on changes in pore connectivity. The link between the two models is not a subject of the work and will be reported in future communications.

Discrete lattice representation of the material microstructure seems to offer the most appropriate modelling strategy for linking the mechanical behaviour to the pore network models. This is a meso-scale approach, where the material is appropriately subdivided into cells and lattice sites are placed at the centres of the cells. From one side, it is a natural solid-phase counterpart to the discrete pore networks. From the other side, discrete lattices allow for studies of distributed damage (microcracks) without constitutive assumptions about crack paths and coalescences that would be needed in a continuum finite element modelling. The deformation of the represented continuum arises from the interactions between the lattice sites. These involve forces resisting relative displacements and moments resisting relative rotations between sites. Two conceptually similar approaches have been proposed to link local interactions to continuum response. In the first one, the local forces are related to the stresses in the continuum cell (see e.g. [5-8]). In the second one, the interactions are represented by structural beam elements (or bonds), the stiffness coefficients of which are determined by equating the strain energy in the discrete and the continuum cell (see e.g. [9-12]). In both cases explicit relations between local and continuum parameters can be established for regular lattices [13]. For irregular lattices, such as those based on random Voronoi tessellation of space, such relations can be established only in an average sense [14]. Regular lattices, however, remain attractive because of the higher computational efficiency and the ability to upscale results.

Most of the previous works on lattice modelling of concrete have used a 2D lattice with hexagonal unit cell [9-12,15]. The reason is that this lattice can be made correspondent to isotropic elastic materials with Poisson's ratio of up to 1/4 in plane strain and up to 1/3 in plane stress [10,14], hence covering a relatively large class of engineering materials. The progress to 3D simulations has been hindered by the fact that the simple 3D lattices cannot be made correspondent to isotropic materials other than materials with zero Poisson's ratio. This is the case for the lattices based on HCP and FCC arrangements

analysed in [14], where each site has 12 coordinating equidistant neighbours. Similar analysis for the simple cubic arrangement shows that again only isotropic elasticity with zero Poisson's ratio can be represented. A bi-regular lattice, based on a truncated octahedron unit cell, has been recently shown to be able to represent any elastic material of practical interest [16]. This quality comes from the presence of bonds with two different lengths, the properties of which can be selected independently, and creates the opportunity for realistic 3D simulations of damage evolution. A pilot study of microcrack evolution in cement with this model and a synthetic distribution of lattice bond strengths has been previously reported [17].

The aim of this work is to demonstrate a relation between a microstructure parameter, the pore size distribution, and a macroscopic behaviour, the stress-strain response of concrete. The success of such a demonstration will support strongly the planned coupling of the mechanical and the transport models. As a first approximation the concrete is considered isotropic elastic and modelled with the novel bi-regular lattice. The concrete non-linear behaviour is thus attributed solely to the emergence of microcracks, represented by failures of lattice bonds. The bond failure is controlled by the forces in the bond and the size of a pore assigned to the bond. The pore sizes are assigned to the bonds according to experimentally determined pore size distribution. The mechanical properties of the concrete under several loading conditions are used to explore the model. The link between the pore size distribution and the stress-strain response of the concrete is clearly demonstrated. The difference between the experimental and the predicted behaviour, particularly at strains close to failure, is attributed to concrete plasticity, currently not included in the model. It is concluded, that for relatively small deformations of the concrete the elastic model presented here will be sufficient for predicting changes in the pore space.

2. Materials and experimental methods

2.1. Pore size distribution

Concrete samples with the composition given in Table 1 were used for this study. Data for the pore size distribution was obtained using X-ray Computed Tomography (X-ray CT). A cylindrical concrete core with a diameter of 28 mm was analysed, using the Nikon 225 kV Custom Bay at the Henry-Moseley X-ray Imaging Facility (HMXIF), The University of Manchester. For data acquisition, a 1.5 mm thick copper filter with accelerating voltage of 160 kV, and a beam current of 200 lA was chosen. 3142 radiographic projections over a 360° rotation were obtained, with an acquisition time of 1000 ms per radiograph. The Metris CT pro software was employed for reconstructing a 3D volume, with corrections for centre shift, beam hardening and noise reduction.

The Aviso Fire segmentation software [18] was used to obtain pore size distribution from the 3D data set. A region of interest (ROI) of 25.2 mm x 18 mm x 18 mm located at the centre of the cylindrical specimen was chosen to reduce possible artifacts that may arise from analysing the free surface and surface breaking pores of the sample. The ROI had a dimension of 1700 x 1200 x 1200 voxels (3D pixels) with a voxel size of ca. 15 im, Fig. 1a. By using the Nyquist-Shannon sampling theorem [19], and assuming at least three pixels per resolvable element, a minimum detectable pore diameter of 31.5 im could be obtained in this study.

The grey-scale value (GSV) distribution of the concrete contained three distinct peaks, with one of those peaks associated with porosity. A minima thresholding routine to differentiate porosity from the cement matrix and the aggregates was used, followed by a sensitivity analysis to inform about the accuracy of the pore threshold, Fig. 1b. Increments of ±25 GSV with respect to the chosen pore threshold produced an error of less than two percent. The segmented voxels corresponding with the porous phase were then analysed using the Avizo XQuant Quantification tool in Avizo Fire. This routine analyses each voxel in the porous phase of the ROI to ascertain any connectivity to other voxels of the same phase. Voxels belonging to the same connected cluster are then grouped together, and the volume of these clusters is used to calculate an equivalent pore size diameter assuming spherical pore morphology.

Fig. 2 shows the results for the probability density, f(c), and cumulative probability, F(x < c), of pore radii. The probability density histogram, Fig. 2a, is determined by an optimised algorithm for bin-size selection [20] yielding a bin size of ca. 2.6 im. The cumulative probability function, Fig. 2b, is obtained with standard median ranking, where for pore radii ordered as c1 6 c2 6 ■■■ 6 cn, the cumulative probability for pores with radii less than ci is given by F(c< ci) = (i - 0.3)/(n + 0.4). The number of pores measured experimentally is n « 8900. The minimum, maximum and average pore radii are depicted in Fig. 2. The cumulative distribution function is used to populate the model described in Section 3 with pore sizes.

2.2. Mechanical tests

All mechanical tests presented in this paper were carried out on cubic concrete specimens, in mac2T, the facility for multi-axial compression of concrete at elevated temperature at The University of Sheffield [21]. The rig comprises three

Table 1

Composition of the concrete used for this study (OPC - Ordinary portland cement, PFA - Pulverized fuel ash).

Constituent OPC PFA Sand Quartz 10 mm Quartz 20 mm Plasticiser SP4 Water

Mass ratio 1 0.33 2.45 1.39 2.78 0.0006 0.56

A.P. Jivkov et al./Engineering Fracture Mechanics xxx (2013) xxx-xxx

Fig. 1. Reconstructed control specimen (a) and segmented porosity within the volume (b).

Fig. 2. Statistical distributions of pore radii in the concrete: (a) probability density; (b) cumulative probability.

independent, orthogonal loading frames with 4MN capacity. The deformations of the specimens are obtained by the means of a system of 6 laser interferometers, measuring the positions of each of the six specimen surfaces. The load is applied via rigid, PTFE coated platens [22]. The concrete was cast in larger 130 mm thick slabs, then cut to 105 mm cubes, and machined (ground) to 100 mm right-regular cubic specimens that were used in the tests. All tests were performed at ambient temperature, using the following loading paths:

Test 1: Unconfined uniaxial compression. This test was performed by increasing one principal stress (r3) monotonically to 80% of peak stress at a rate of 14.7 MPa/min, followed by strain control loading at 60 im/min to the post-peak region. The other two stresses were kept constant at 1.0 MPa in order to maintain contact between the loading platens and the specimen, needed for the deformation measurements. Recorded stress-strain response is illustrated in Fig. 3a. Test 2: Extension meridian loading under confinement. This test was performed in two stages: (i) loading hydrostatically to confinement p = r1 = r2 = r3 = 61 MPa, (ii) keeping the minor stress r2 constant at 61 MPa (confinement) and

Fig. 3. Experimental stress-strain behaviour of concrete: (a) unconfined uniaxial compression; (b) extension meridian loading under confinement r2 = 61 MPa (Test 2); (c) shear meridian loading under confinement r2 = 61 MPa (Test 3).

A.P.Jivkov et al./Engineering Fracture Mechanics xxx (2013) xxx-xxx

increasing the two major stresses o^ = a3 = p + Da to the peak, at the same stress/strain rates as in Test 1. Recorded stress-strain response is illustrated in Fig. 3b.

Test 3: Shear meridian loading under confinement. This test was performed by following the same procedure as that in Test 2, except that in the second stage the major and the intermediate stress were different: a3 = p + Da and a1 = p + Da2. Recorded stress-strain response is illustrated in Fig. 3c.

Test 4: Cyclic unconfined uniaxial compression. This test was performed by following the same procedure as that in Test 1, except that the major stress (a3) was applied in three loading-unloading cycles, between 1 MPa and 18, 30 and 42 MPa, before loading to the peak (at 60 MPa).

Test 5: Hydrostatic compression. This test was performed in three cycles, by loading a1 = a2 = a3 between 1 and 110 MPa, 220 MPa and 330 MPa.

The elastic properties of the material were calculated from the unloading branches of cyclic uniaxial (unconfined) compression tests, (Test 4: Modulus of elasticity E = 46 GPa, and Poisson's ratio v = 0.27), and hydrostatic compression tests (Test 5: Bulk modulus K = 33 GPa).

3. Modelling and simulations

3.1. Lattice model

Fig. 4 illustrates the lattice model used in this work. The unit cell, shown in Fig. 4a is a truncated octahedron - a solid with six square and eight regular hexagonal boundaries. The 3D space can be compactly tessellated using such cells, with each cell representing a material meso-scale feature, e.g. grain, in an average sense. This representation is supported by physical and statistical arguments [16]. A discrete lattice is formed by placing sites at the centres of the cells and connecting each site to its 14 nearest neighbours; example is shown in Fig. 4b. The lattice contains two types of bonds. Bonds denoted by B1 in Fig. 4 are normal to the square boundaries of the unit cell and form orthogonal set. For convenience this set is made coincident with the global coordinate system and B1 are referred to as principal bonds. Bonds denoted by B2 in Fig. 4 are normal to the hexagonal boundaries of the unit cell. The hexagons lie on the octahedral planes with respect to the selected system. Hence B2 are referred to as octahedral bonds. If the spacing between sites in the three principal directions is denoted by L, bonds B1 have length L1 = L, and bonds B2 have length L2 = p3L/2.

One important point about the proposed lattice is that it does not permit a simple closed form solution for the relation between bond properties and continuum elastic constants. In the general case, the behaviour of the two distinct types of bonds, B1 and B2, needs to be represented by eight parameters - normal, shear, twisting and bending stiffness for each bond type. With analysis based on equivalence of the energy in the continuum and the discrete unit cell, e.g. following [13], it can be shown that the proposed lattice represents a micropolar material with cubic elasticity. This is not surprising as the unit cell of the proposed model is the Voronoi cell (or the first Brillouin zone) of the face-centred cubic crystals. With analysis of a homogeneous displacement field, i.e. with no relative rotations between sites, the three constants of the cubic elasticity can be related to the four linear stiffness coefficients (normal and shear) of the bonds. Hence the model is over-determined, i.e. a particular set of continuum elastic constants can be achieved with infinite number of sets of bond stiffness coefficients. With analyses of more complex deformation modes involving relative rotations between sites, relations between the linear and the twisting/bending stiffness coefficients can be obtained. If such analyses are based on deformation energy functional dependent on strains only (classical continuum mechanics used in [13]) the problem remains over-determined. One possible explanation is that the proposed lattice arrangement is genuinely "micropolar", i.e. the stiffness coefficients of the bonds

Fig. 4. Lattice illustration: (a) unit cell (truncated octahedron) showing the site with 14 coordinating bonds - six in principal directions normal to squares (B,) and eight in octahedral directions normal to hexagons (B2); (b) discrete lattice of beam elements.

cannot be uniquely determined from classical strain energy potential. This is a subject of ongoing work and will be reported in a future publication.

The lattice, however, has been demonstrated to produce any predefined elastic response, including isotropic elasticity for a large class of materials [16]. This demonstration uses a "global" approach in the sense that the strain energy stored in the entire lattice is compared to the strain energy stored in a continuum with corresponding volume. It is assumed that for sufficiently large lattices the global approach provides satisfactory approximation with intrinsic micropolarity averaged over the volume. With this assumption it is convenient (and sufficient) to represent the bonds of the lattice with structural beam elements of circular cross sections, with R1 and R2 denoting the radii of beams B1 and B2, respectively. The beams are clamped at the lattice sites, i.e. they transfer three forces and three moments representing the resistance of the material to relative displacements and rotations between adjacent solid cells or grains. Assuming local isotropy, the modulus of elasticity, Eb, and Poisson's ratio, vb, of the two types of beams should be the same. The four parameters, R1/L, R2/L, Eb, and vb can be calibrated so that the lattice produces the required isotropic elastic response. For the concrete studied in this work with E = 46 GPa and v = 0.27, the calibration following the procedure in [16] yields R1/L = 0.2; R2/L = 0.32; Eb = 90 GPa; and vb = 0.4. These parameters are used for the simulations reported in the Section 4. The commercial software Abaqus [23] with Euler-Bernoulli beam formulation has been used. For the calibrated radius to length ratio of the beams, the use of the Timoshenko beam formulation would not offer noticeable improvement because the benefits are typically for ratios up to 1/8.

Pores of different sizes from measured X-ray CT data pore size distribution are assigned to the lattice bonds. A generator of uniformly distributed random numbers 0 6 r <1 is used to assign experimentally determined pore sizes to individual bonds. For each bond a random number r is generated and the assigned pore radius is calculated from c = F_1(r), where F(c) is the cumulative probability of pore sizes given in Fig. 2b. This is a standard statistical process that ensures that the distribution of pore sizes in the model comes from the same population as the experimentally determined pore sizes. Fig. 5 illustrates a fragment of the model with distributed pores. The cellular structure is shown in order to introduce the length relative to which the pore sizes are to scale. With respect to the cellular structure the pores reside at cell boundaries, i.e. interfaces between grains. The bonds of the corresponding lattice model are also depicted (not to scale) in order to show that pores reside at bond centres.

3.2. Failure criterion

Damage in the lattice model is introduced by removal of bonds. The criterion for bond failure is based on the forces and moments in the beam elements. This is formulated by

Nf + Sf + Tf + Mf P (>

where N and S are the normal and shear forces in the beam, respectively; T and M are the twisting and the bending moments, respectively; Nf, Sf, Tf, and Mf are critical values of these forces and moments. The normal force is taken with its sign - positive for tension and negative for compression. The shear force and the bending moment in Eq. (1) are obtained from the values in the two directions normal to the beam axis using the square root of squares rule. This criterion has been suggested previously to account for the contribution of all deformation modes to failure [24]. Thus, a bond is permitted to fail under pure extension, pure shear, pure twist and pure bending separately, as well as under the combined action of the forces and moments. Taking only the first and the fourth term in Eq. (1) is principally equivalent to some of the previously used criteria, e.g. [9], where shear failure is not accounted for. Introducing the second and the third term allows for shear failure similarly to

Fig. 5. Segment of model illustrating pores of different sizes distributed to cell boundaries and corresponding lattice bonds. Pore sizes shown are to scale with the cell size and follow the experimentally determined pore size distribution. Lattice bond diameters are not to scale.

the way used in [11]. Note that in the current model failure may occur in a compressed bond provided that the other forces are sufficiently large to fulfil the criterion. Note also that in a bond with tensile normal force, the contribution of the shear force and the moments can lead to failure of the bond before the normal force reaches its critical value. This provides an interaction between the different forces that allows for pore failure under the combined action of normal and shear stresses. If the sum of the first and fourth terms in Eq. (1) is larger than the sum of the second and third term, the failure is dominated by the normal stresses and we shall call this failure by separation. If the inverse is true, the failure is dominated by the shear stresses and we shall call this failure by sliding.

The failure parameters for the four modes, Nf, Sf, Tf, and Mf, can be related with the following argument [24]. For a beam of circular cross section of radius R, the failure stress corresponding to Nf is Of = Nf/(p R2); a tensile failure stress. The maximum stress due to bending is amax = 4 M/(p R3). This equals the tensile failure stress for Mf = Nf R/4. Similarly, the failure stress corresponding to Sf is Sf = Sf/(p R2); a shear failure stress. The maximum stress due to torsion is smax = 2 T/(p R3) and this equals the shear failure stress for Tf = Sf R/2. Thus the failure criterion given by Eq. (2) requires two material parameters - the tensile and shear failure stresses, or equivalently failure forces. For convenience we introduce the parameter g = Sf/Nf = Sf/of. Typical values for the shear to tensile strength ratio g are between one and two for quasi-brittle materials, with larger values yielding more brittle behaviour in uniaxial compression tests. This will be demonstrated by a parametric study in Section 4.

Key feature of the failure model is the relation between the tensile failure strength of a bond, Of, and the size of the pore assigned to the bond. This is based on the average stress criterion, proposed originally for circular holes (and cylindrical voids in 3D) [25]. This criterion accounts for the size effect on the failure strength by recognizing that the difference between two voids of identical shapes but different sizes is in the stress profiles ahead of the void surfaces. The average stress is an integral of the stress profile from the void surface to a material dependent distance, d, divided by this distance. The criterion states that failure occurs when the average stress reaches a critical value, a0. For a given material distance, the average stress is larger for a larger pore, providing the required differentiation between pore sizes where a larger pore fails at lower applied stress. The material dependent parameter, d, is a function of, for example, pore space, aggregate size, cement grains, or hydrated/un-hydrated regions present in the cement matrix. The critical strength, a0, can be interpreted as the "ideal" or "theoretical" tensile strength of the material, i.e. the strength without defects present. For convenience we introduce the non-dimensional parameter n = c/(c + d), where c is the pore radius. In terms of n the tensile failure strength of a bond with a spherical void is related to the "ideal" strength via

Of =_2__(2)

a0 2 + n(1 + n)[2+28(1 + n2)] w

Eq. (2) is found by integrating the stress profile away from the surface of a spherical void over distance d and dividing by this distance [26]. Thus the failure model is based on three material dependent parameters: the material length scale, d; the "ideal" tensile strength, a0; and the shear to tensile strength ratio, g. The material length scale is expected to be of the order of the smallest pore size. Parametric study, shown in Section 4, suggests that length scales around the smallest pore size have effect on the macroscopic response up to roughly the average pore size after which the effect is negligible.

3.3. Simulation details

With respect to a coordinate system (X1,X2,X3), coincident with the principal bonds, a model of a cubic region (20L,20L,20L) has been used. The corresponding lattice contains 17261 sites and 113260 bonds - 49260 of type B1 and 64,000 of type B2. Note that for this lattice, each of the boundary planes X1 =0, X1 = 20L, X2 = 0, X2 = 20L, X3 = 0, X3 = 20L contains 21 x 21 sites. Three loading cases have been simulated: the unconfined uniaxial compression (UC); the extension meridian loading (EM); and the shear meridian loading (SM). The lattice has been loaded via prescribed displacements corresponding to the strains measured experimentally. Let (U1, U2, U3) are the displacements and (F1,F2,F3) are the reaction forces of sites on the lattice boundaries with respect to the coordinate system (X1,X2,X3). For all loading cases the following boundary conditions have been fixed: U1 = 0 for sites on X1 = 0; U2 = 0 for sites on X2 = 0; U3 = 0 for sites on X3 = 0.

For UC, the conditions for the remaining boundaries are: F1 = 0 (free sites) on X1 = 20L; F2 = 0 (free sites) on X2 = 20L; U3 = 20Le3(t) for sites onX3 = 20L, where e3(t) is the experimentally measured compressive strain evolution. The macroscopic stress in the loading direction is resolved as the ratio between the total reaction force at sites on X3 = 20L and the boundary area, i.e. a3 = EF3/400L2. The macroscopic strains in the lateral directions are taken as the ratios between the average displacements of sites on X1 = 20L on X2 = 20L, and the model length, i.e. e1 = SU1/(212 x 20L) and e2 = EU2/(212 x 20L).

For EM and SM, the lattice has been initially subjected to hydrostatic compression with applied displacements: U1 = 20Le(t) for sites on X1 = 20L; U2 = 20Le(t) for sites on X2 = 20L; U3 = 20Le(t) for sites on X3 = 20L, where the compressive strain e(t) has been increased until the macroscopic boundary stresses a1 = a2 = a3 reached the confining pressure p = 61 MPa as in experiment. The reaction forces at sites on X2 = 20L are recorded for the next loading step as F2(p). Within this, the boundary conditions have been changed to: U1 = 20Le1(t) for sites on X1 = 20L; U3 = 20Le3(t) for sites on X3 = 20L; F2 = F2(p) for sites onX2 = 20L. Here, e1(t) and e3(t) are the experimentally measured compressive strain evolutions. The macroscopic stresses a1 and a3 are determined as above, while a2 = 61 MPa. The macroscopic strain e2 is determined as above.

The simulations have been performed with an in-house code controlling the failure of bonds and interfacing with Abaqus, which is used for solving the lattice equilibrium after each load increment. Variable load increments have been used, so that a single failure event occurs within one increment. This is based on the standard cut-back algorithm used in plasticity. Together with the emergent macroscopic behaviour (stress versus strain), the failure events and their nature (separation versus sliding) have also been recorded for the bonds in all lattice directions. The model has been calibrated using the unconfined uniaxial compression loading case. We define the "ideal" tensile strength, r0, to be the value for which the maximum stress in the simulations equals the maximum stress in the experiment, see Fig. 3a. This depends on the selected material length scale, d, and the shear to tensile strength ratio, g, which is shown by parametric studies in Section 4. A particular selection of the three parameters is used for the simulations of the extension and the shear meridian loading cases. This is given by d = 40 im and g = 2, for which the calibrated "ideal" tensile strength is r0 = 28 MPa. The selection is related to the outcomes from the parametric studies under unconfined uniaxial compression.

4. Results

All results presented in this section are obtained with one and the same distribution of pore sizes within the bonds. Statistical analysis of the model responses for various distributions is outside the scope of the present work. In all figures presented, the strains and stresses are assumed positive if compressive and negative if tensile.

4.1. Unconfined uniaxial compression: effect of material length scale

Fig. 6 shows the results of the simulations performed to study the effect of the material length scale, d, with fixed shear to tensile strength ratio and ideal tensile strength as depicted. The macroscopic stress in the direction of the applied displacement is plotted versus the macroscopic strain. The experimental stress-strain behaviour is shown for comparison with broken line. The four values of d, depicted in the figure, are selected from less than the minimum pore size (15 im) up to the average pore size (50 im). The simulated behaviour follows closely the experiment for sufficiently small strains. The effect of smaller length scales is in the reduction of the predicted peak stress after which softening occurs. For the results in Fig. 6, the ideal strength is selected such that the peak stress for the case d = 40 im equals the experimental one. Hence, a smaller length scale requires a larger ideal tensile strength value to match the experimental behaviour. However, the effect of the length scale diminishes for d >50 im (the average pore size) as seen by comparing the curves for d = 40 im and d = 50 im.

The microcracking within the volume studied is initially randomly distributed. It is therefore convenient to represent the development of damage as fractions of failed bonds relative to the total number of bonds in the model. Fig. 7 shows the fractions of bonds failed by separation (a) and by sliding (b) for the four length scale cases presented in Fig. 6. The point of peak stress recorded for each length scale case is shown with a star on the corresponding curve. These curves show that the evolution of damage commences with rapid increase of separation cracks which flattens out with the onset of sliding cracks; compare Fig. 7a with Fig. 7b. The number of separation cracks at peak stress (onset of softening behaviour) is nearly independent of the material length scale selected, Fig. 7a; the peak stress occurs at the inflection points of the curves. The post-peak (softening) behaviour is characterised with a rapid increase of sliding microcracks, Fig. 7b, and continuing but moderate increase of separation microcracks, Fig. 7a. There is a small effect of the material length scale on the sliding failures at peak stress, Fig. 7b, however, this diminishes with increasing the length scale. Although the post-peak behaviour is still affected by the material length scale, we choose the case d = 40 im as basic for the subsequent parametric studies and analysis.

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 ¿3XIO-3

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5

¿3X10-3 (b)

Fig. 7. Evolution of damage (microcracks) under unconfined uniaxial compression, represented with fractions of failed bonds out of all bonds in the model: (a) separation microcracks and (b) sliding microcracks. The points of recorded peak stress are shown with stars on the four length scale cases curves.

To further illustrate the nature and orientation of the distributed microcracks, the separation and sliding failures normal to the principal and octahedral directions are considered for the basic case with d = 40 im, g = 2, and r0 = 28 MPa. Fig. 8 shows the results for separation microcracks (a) and sliding microcracks (b). The principal bonds in the unconfined directions, i.e. those parallel to e1 and e2, and denoted by B1(e1) and B-i(e2), respectively, develop exclusively separation microcracks and fail rapidly during the initial stages of straining; see Fig. 8a and b. These are followed by rapid growth of separation and sliding failures on the octahedral planes, denoted by B2. A small number of separation and sliding microcracks develop also in the principal bonds parallel to the loading direction, denoted by B1(e3). These form around and after the peak during the softening and are apparently due to excessive bending contribution in the failure criterion.

Fig. 9 shows the prediction of the model for both the compressive and the lateral strains compared to the experimental measurement shown in Fig. 3a. The development of the strain in the unconfined direction is in a very good agreement with experiment up to strains very close to the peak of the response. This supports strongly the capability of the proposed model to simulate the pre-peak macroscopic behaviour based on the elastic calibration of the lattice and the development of distributed failures governed by the sizes of the pores. The post-peak behaviour of the compressive and the lateral strains cannot be reproduced correctly. This suggests that other types of non-linearities become dominating, e.g. plastic behaviour and geometry effects during testing.

4.2. Unconfined uniaxial compression: effect of shear to tensile strength ratio

Fig. 10 shows the results of the simulations performed to study the effect of the shear to tensile strength ratio, g, with fixed material length scale and ideal tensile strength as depicted. Typical values of g are studied, including g = 2 shown before in Fig. 6. The effect of g is similar in nature to the effect of material length scale, see Fig. 5. The peak stress is reached at lower strains for lower shear strength values. This suggests more brittle behaviour for materials with lower shear to tensile strength ratio. For known g <2, the model can be recalibrated to find the ideal tensile strength, r0, for which the predicted peak stress equals the experimental one.

A.P. Jivkov et al./Engineering Fracture Mechanics xxx (2013) xxx-xxx

ra 1E-4

0.0 0.5

: B^) ^__

-- Peak stress

2.0 £3 x 10 3

= 0.01

No sliding in B, normal to e1 and e2

Peak stress --

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 Î3X10"3

Fig. 8. Evolution of microcracks under unconfined uniaxial compression, represented with fractions of failed bonds out of all bonds in the model: (a) separation microcracks on principal and octahedral planes and (b) sliding microcracks on principal and octahedral planes.

Fig. 9. Simulated stress-strain behaviour of concrete under unconfined uniaxial compression: strains in compression and unconfined directions. Results for the basic selection of model parameters depicted.

The development of microcracks for different values of g is illustrated in Fig. 11. The curves in Fig. 11a shows the fraction of bonds failed by separation and the curves in Fig. 11b shows the fraction of bonds failed by sliding. These show that damage starts with rapid increase of separation failures, similarly to the observation in Fig. 7, which slows down after the onset of sliding failures. The increase of sliding failures dominates the response up to the peak stress with very little contribution of separation failures. The post-peak (softening) response is characterised by continuing rapid increase in sliding and moderate

\ \ -------- 7 = 2

/ \ V \ /7=1.75

/ '/= 1 N

71= 1.5

7= 1.25

d= 40 |jm; a0 = 28 MPa

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0

i3x10-3

Fig. 10. Simulated stress-strain behaviour of concrete under unconfined uniaxial compression: effect of shear to tensile strength ratio (depicted on curves) for fixed material length scale and ideal tensile strength (depicted on graph).

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 ^xlO"3

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5

£3x103 (b)

Fig. 11. Evolution of microcracks under unconfined uniaxial compression, represented with fractions of failed bonds out of all bonds in the model: (a) separation microcracks and (b) sliding microcracks. The points of recorded peak stress are shown with stars on the four length scale cases curves.

increase in separation failures. Notably, a reduction of the shear strength, i.e. decreasing g, leads to decreasing contribution of separation failures and increasing contribution of sliding failures near the peak stress and during softening. For sufficiently small strains, e.g. e3 < 1 x 10~3, the effect of g is negligible, Fig. 11a. The behaviour is dominated by separation failures normal to the unconfined planes of the model and on octahedral planes, see Fig. 8a.

4.3. Complex loading: extension and shear meridians

Fig. 12 shows the simulated response of the concrete under the extension meridian loading, compared to the experimentally measured response. The initial confining pressure and the model parameters used are depicted. The predicted responses

e2 (mod) e3 = e, (mod)

£2 (exp) \ L // h = «I (exP)

\f cf = 40 |jm

- p = 61 MPa ! 77=2

<j0 = 28 MPa

-2-10123

Strain x 10"2

Fig. 12. Simulated stress-strain behaviour of concrete under extension meridian loading: strains in compression and unconfined directions. Experimental curves (extension meridian loading under 61 MPa confinement) included for comparison. Model parameters and initial confining pressure depicted.

0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 s3x10-2

0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25

Fig. 13. Evolution of microcracks under extension meridian loading, represented with fractions of failed bonds out of all bonds in the model: (a) separation microcracks on principal and octahedral planes and (b) sliding microcracks on principal and octahedral planes.

A.P.Jivkov et al./Engineering Fracture Mechanics xxx (2013) xxx-xxx

are in very good agreement with experiment for sufficiently small, but well above the confining, stresses. However, for larger strains the simulations show "stiffer" and more brittle response. This is an indication that the non-linearity introduced by the elastic microcracking model is not sufficient for describing the macroscopic behaviour at high stresses and strains. However, a significant part of the non-linear response can be attributed to microcracking.

Fig. 14. Simulated stress-strain behaviour of concrete under shear meridian loading: strains in all principal directions. Experimental curves (shear meridian loading under 61 MPa confinement) included for comparison. Model parameters and initial confining pressure depicted.

Fig. 15. Evolution of microcracks under shear meridian loading, represented with fractions of failed bonds out of all bonds in the model: (a) separation microcracks on principal and octahedral planes and (b) sliding microcracks on principal and octahedral planes.

Fig. 13 shows the evolution of microcracking under the extension meridian loading leading to the emergent stress-strain response shown in Fig. 12. The notations in Fig. 13 are identical to the notations in Fig. 8. Notably, damage commences with rapid development of separation, Fig. 13a, and sliding, Fig. 13b, microcracks in bonds in the direction of constant stress, B1(e2). The damage continues with increase of the separation failures in the same bonds, Fig. 13a, followed by rapid increase of sliding failures on octahedral planes, B2, Fig. 13b. This leads to the peak in the simulated stress response, see Fig. 12, after which the softening is characterised by increases in separation and sliding failures in the bonds in the directions of continuous loading, B1(e1) and B1(e3), and on octahedral planes, B2. Although the simulations may not be representative for higher strains, the results in Fig. 13 provide useful information on the nature of microcracking. The observed behaviour is physically realistic and illustrates the potential of the model.

Fig. 14 presents the simulated response of the concrete under the shear meridian loading, compared to the experimentally measured response. The agreement with experiment for this loading case is somewhat better than the agreement for the extension meridian, Fig. 12. In particular, the strain in the direction of constant stress, e2, is in excellent agreement up to very high stresses. The responses are less stiff and brittle for all principal strains, suggesting that the contribution of the elastic microcracking to the overall non-linear response is larger under shear meridian loading than under extension meridian loading. The contribution of plasticity remains important at increasing strains.

Fig. 15 illustrates the evolution of microcracking under the shear meridian loading leading to the emergent stress-strain response shown in Fig. 14. Damage commences with rapid development of separation, Fig. 15a, and sliding, Fig. 15b, microcracks in bonds in the direction of constant stress, B1(e2). This is followed by rapid increases in sliding failures on octahedral planes, B2, and in bonds in the direction of the smaller increasing strain, B1(e1), Fig. 15b. In parallel, the separation failures in B1(e2) continue to grow together with emerging separation failures in B1(e1), Fig. 15a. This leads to the peak in the simulated stress response, see Fig. 14, after which the softening is characterised by increases in separation and sliding failures in all bonds, including those in the direction of the larger increasing strain, B1(e3). Comparison between Fig. 15 and Fig. 13, shows that the shear meridian loading produces a more complex microcrack population. The peak stress is reached due to separation and sliding failures in two principal directions in addition to sliding failures on octahedral planes.

The biggest differences between the complex loading cases and the unconstrained uniaxial compression are: (i) for the confined cases separation failures on octahedral planes do not occur prior to softening (see B2 in Figs. 13a and 15a), while for the unconfined case such failures contribute substantially to peak stress (see B2 in Fig. 8a) and (ii) for the unconfined case sliding failures on the planes parallel to the load do not occur (see B1(e1) and B1(e2) in Fig. 8b), while for the confined cases such failures contribute to the peak stress (see B1(e1) and B1(e2) in Figs. 13b and 15b).

5. Discussion

The model presented in this work is based on a linear elastic behaviour at the material meso-scale, i.e. linear elastic bonds. The non-linearity in the simulated macroscopic responses is solely due to elimination of bonds representing the emergence of microcracks. The criterion for bond elimination (failure) is based on the local forces and moments and on the size of a pore allocated to the bond. The pore size effect is accounted for using the average stress over a material distance, with failure occurring when this stress reaches a critical "ideal" strength.

The material distance, or length scale, is not known. However, it has been demonstrated that the effect of this parameter on the predicted response is significant only for values of the order of the smallest pores and diminishes quickly for values above the average pore size, see Fig. 6. The results suggest that smaller material length scales yield more brittle macroscopic responses, as the peak stress (onset of softening) is reached at lower strains. The pore size distribution measured experimentally, Fig. 2, represents the tail of the real distribution as pores of sizes below 15.75 im have not been resolved. However, the average pore size will not change noticeably if the missing head of the distribution is included, e.g. by assuming a Weibull-continuation of the measured data. It could be reduced from the current 50 im down to at most 40 im. Hence, the results presented in this work will not change substantially with the inclusion of smaller pore sizes, provided that the material length scale is larger than the average pore size. However, if the material length scale is known from experiment to be smaller than the average pore size, the model can be recalibrated to find the "ideal" tensile strength for which the simulated response matches the experimental one.

The other parameter controlling the macroscopic behaviour is the shear to tensile strength ratio. Experimental value for this parameter is not available for the concrete studied in the work. Parametric studies with the typical values for quasi-brittle materials performed illustrate the effect of the shear strength. It has been shown that smaller shear strength values produce more brittle macroscopic responses, Fig. 10. Notably, for a known shear to tensile strength ratio, the model can be recalibrated to find the corresponding "ideal" tensile strength.

The "ideal" tensile strength is a model parameter representing the uniform tensile stress required to break a material volume with no defects present. Therefore it cannot be related to the experimentally determined tensile strength of the concrete in a simple manner. As the results suggest, it depends on the material length scale and the shear strength, Figs. 6 and 10. The value used in the work, r0 = 28 MPa, can be considered as a lower bound for the parametric studies performed. This is because any reduction of the shear strength, as well as any reduction of the material length scale, would require an increase of r0 in order to simulate the experimental macroscopic response. Therefore the "ideal" tensile strength remains a calibrating parameter at the meso-scale. The model presented considers homogenised meso-scale properties, i.e. the solid

phases of the concrete and their interactions are not accounted for separately. In principle, the model can be extended to include different r0 for the interactions between different phases, e.g. by assigning appropriate values to different bonds. However, for small differences in r0 for different phases, the effect will be negligible when compared to the effect of the sizes of the distributed pores. Nevertheless, this possibility should be explored in the future, in particular if the sizes of the pores are found to be correlated to the solid phases, e.g. part of the distribution is found predominantly in one phase.

The model realism is strongly supported by the results presented for the nature and orientation of microcracks. In the case of unconfined uniaxial compression, Figs. 7, 8 and 11 demonstrate that damage initiates by separation failures normal to the unconfined directions followed by separation failures on octahedral planes, Fig. 8a. The rapid increase in separation microcracks is interrupted by the occurrence of sliding failures on octahedral planes. The onset of sliding failures can be considered as the compressive strength of the concrete, because the fast increase of these leads quickly to the peak stress in the simulated response. The post-peak behaviour is a well-defined material softening due to continuing rapid increase in sliding failures and a more moderate increase in separation failures. This observations correlate very well with the expected development of microcracking. The resulting macroscopic stress-strain response in the loading and lateral directions, Fig. 9, is in excellent agreement with the experimental measurement for strains up to the onset of sliding failures. This is despite the use of linear elastic behaviour at the meso-scale. This suggests that the unconfined compression behaviour is controlled predominantly by microcracking in the region before the compressive strength is reached. Within the peak stress and post-peak regions, characterised by rapid growth of sliding failures in the elastic model, the predictions for loading and lateral strains deviate from the experiment. This is an indication that the measured response cannot be attributed to brittle microcracking alone; it becomes dominated by other factors, such as plasticity and creep at the meso-scale.

Our conclusions for the model realism are qualitative only, due to lack of experimental evidence for the microcrack formation prior to peak load. A rare experimental study that reports direct observations of damage in concrete using X-ray mic-rotomography [27], presents only macro-cracks at peak stress, under cycling uniaxial compression. Attempts have been made to detect localization zones by measuring velocities of ultrasonic waves [28], monitoring the displacement field by Electronic Speckle Pattern Interferometer [29], and using acoustic emission [30,31]. In principle we could make use of these techniques in the future to verify more quantitatively the uniaxial compression results. However, the techniques are either impossible or very difficult to implement in multiaxial compression tests, where all surfaces of the specimens are covered by the loading platens.

The results for the complex loading cases, extension meridian and shear meridian after hydrostatic compression, provide further evidence that the model generates qualitatively physically realistic evolution of microcrack populations. The predicted stress-strain responses are also in good agreement with experiment at lower stresses. For higher stresses it is clear that the combination of elastic bond and brittle damage cannot account for stress strain behaviour under confinement. In this regime the local behaviour is predominantly plastic and potentially time-dependent. There are no models published in literature that can successfully simulate post peak behaviour recorded in experiments carried out in multiaxial compression. Moreover, experiments in the post-peak region show that the test results are sensitive to specimen boundary conditions and size of specimens [32-34], as well as applied stress history [35,36]. Despite improvements in boundary conditions, the softening curve still largely depends on size and shape of specimens and their interaction with the loading platens [37].

The differences between the simulated and the experimental stress-strain behaviours confirm the limitations imposed by the simplifications of elastic-brittle bond modelling. The only models that show close agreement of pre-peak stress-strain behaviour under multiaxial stress states are those based on plasticity. Plasticity models based on the Willam-Warnke failure criterion, which include hardening-softening laws, a non-associated flow rule and modifications of the original failure envelope, e.g. [38-43], do provide good agreement with macroscopically observed stress-strain behaviour and good predictions of peak stress. These models, however, depend on non-physical parameters that do not take into account changes in the fabric of the material, and hence are not suitable for coupling with transport equations. Even the continuum damage models, e.g. [44,45], or models based on combinations of plasticity and damage, e.g. [46-48], treat damage as a set of variables characterising the material degradation, which are not explicitly related to the actual changes of the microstructure.

Generally the plastic deformations would reduce the forces and moments in bonds at a given macroscopic strain, which will lead to more moderate development of microcracks and thence to a predicted response closer to the experiment. However, the results illustrate the important contribution of microcracking to the non-linear response of the concrete. This can be used in developing a plasticity model at the meso-scale by comparing the elastic to the elastic-plastic microcracking evolution. Further, the complex loading tests are performed to stress levels significantly higher to those expected in concretes in an underground repository. After sealing such a repository, the principal stresses in the cement-based components will be comparable to the stresses in the near field geology. For a repository at a depth of 500-1000 m, for example, the maximum principal stress is of the order of 15-30 MPa; the other principal stresses being smaller. Hence the concretes will be subjected to stress triaxiality within the elastic range. This could be amended by internally produced stresses due to gas generation for example and externally produced stresses due to small land motions. Unless a catastrophic event, such as an earthquake, occurs, the development of stresses would remain in the elastic region. For this case the results of this work are sufficiently representative. This suggests that the elastic microcracking model can be used to inform pore space models about changes in connectivity leading to predictions for evolution of transport properties. The current model, however, does not include any time-dependent effects on the material properties at the meso-length scale. The incorporation of such effects in the bond properties is in our future plans as this is essential for the long time-scales of repository operation. We intend to

maintain the elastic-brittle behaviour of bonds, considering the relatively low stresses involved, but alter the elastic constants and the failure criterion as a result of time-dependent chemical changes. Within the current framework this means changes to the elastic modulus of bonds and to the critical material length-scale controlling failure around pores.

6. Conclusions

• 3D lattice modelling is a promising approach for correlating microstructural properties to the macroscopic behaviour of concretes.

• Clear link between concrete porosity, in terms of experimentally determined pore size distribution, and the emergent stress-strain response is demonstrated.

• Non-linear response prior to the compressive strength of the concrete is dominated by brittle microcracking formed by local material separation.

• Other material non-linearity, e.g. plasticity and creep, governs the behaviour after the onset of sliding microcracks in the current elastic model, for both uniaxial and triaxial loading cases.

• The elastic model is sufficient to inform pore network models about changes in pore connectivity for the stress levels expected in concrete in underground repository.

• Further development is needed to incorporate time-dependent effects in bond behaviour in order to investigate the concrete behaviour at the time-scales of repository operation.

Acknowledgements

A.P.J. and D.L.E. are supported by BNFL (UK) endowment for research on radioactive waste disposal. R.S. is funded with an EPSRC studentship through the Nuclear FiRST DTC Grant EP/G037140/1. Both supports are highly appreciated. The authors also express their gratitude to HMXIF (The University of Manchester) for the use of X-ray tomography equipment. Funding is acknowledged from the EPSRC for the X-ray imaging facility under EP/F007906/1 and EP/F028431/1.

References

[1] Bernabé Y, Li M, Maineult A. Permeability and pore connectivity: a new model based on network simulations. J Geophys Res 2010;115:B10203.

[2] Ioannidis MA, Chatzis I. Network modelling of pore structure and transport properties of porous media. Chem Engng Sci 1993;48:951-72.

[3] Blunt MJ, Jackson MD, Piri M, Valvatne PH. Detailed physics, predictive capabilities and macroscopic consequences for pore-network models of multiphase flow. Adv Water Resour 2002;25:1069-89.

[4] Jivkov AP, Hollis C, Etiese F, Withers PJ. A novel architecture for pore network modelling with applications to permeability of porous media. J Hydrol 2013;486:246-58.

[5] Bazant ZP, Oh BH. Microplane model for progressive fracture of concrete and rock. J Engng Mech 1985;111:559-82.

[6] Carol I, Jirasek M, Bazant Z. A thermodynamically consistent approach to microplane theory. Part I. Free energy and consistent microplane stress. Int J Solids Struct 2001;38:2921-31.

[7] Cusatis G, Bazant ZP, Cedolin L. Confinement-shear lattice CSL model for fracture propagation in concrete. Comput Meth Appl Mech Engng 2006;195:7154-71.

[8] Cusatis G, Pelessone D, Mencarelli A. Lattice discrete particle model (LDPM) for failure behaviour of concrete. Part I: Theory. Cem Concr Compos 2011;33:881-90.

[9] Schlangen E, van Mier JGM. Experimental and numerical analysis of micromechanisms of fracture of cement-based composites. Cem Concr Compos 1992;14:105-18.

[10] Griffiths DV, Mustoe GGW. Modelling of elastic continua using a grillage of structural elements based on discrete element concepts. Int J Numer Methods Engng 2001;50:1759-75.

[11] Chang CS, Wang TK, Sluys LJ, Van Mier JGM. Fracture modeling using a micro structural mechanics approach - Part I. Theory and formulation. Engng Fract Mech 2002;69:1941-58.

[12] Karihaloo BL, Shao PF, Xiao QZ. Lattice modelling of the failure of particle composites. Engng Fract Mech 2003;70:2385-406.

[13] Wang Y, Mora P. Macroscopic elastic properties of regular lattices. J Mech Phys Solids 2008;56:3459-74.

[14] Grassl P, Jirasek M. Meso-scale approach to modelling the fracture process zone in concrete subjected to uniaxial extension. Int J Solids Struct 2010;47:957-68.

[15] Liu JX, Zhao ZY, Deng SC, Liang NG. Numerical investigation of crack growth in concrete subjected to compression by the generalised beam lattice model. Comput Mech 2009;43:277-95.

[16] Jivkov AP, Yates JR. Elastic behaviour of a regular lattice for meso-scale modelling of solids. Int J Solids Struct 2012;49:3089-99.

[17] Jivkov AP, Gunther M, Travis KP. Site-bond modelling of porous quasi-brittle media. Miner Mag 2012;76(8):94-9.

[18] AVIZO Software - <http://www.vsg3d.com/avizo/fire>.

[19] Jerri AJ. The shannon sampling theorem - its various extensions and applications: a tutorial review. Proc IEEE 1977;65:1565-95.

[20] Shimazaki H, Shinomoto S. A method for selecting the bin size of a time histogram. Neural Comput 2007;19:1503-27.

[21] Petkovski M, Crouch RS, Waldron P. Apparatus for testing concrete under multiaxial compression at elevated temperature (mac2T). Exp Mech 2006;46:387-98.

[22] Petkovski M. Effects of stress during heating on strength and stiffness of concrete at elevated temperature. Cem Concr Res 2010;40:1744-55.

[23] ABAQUS. User's Manual, Version 10.1 ed., Dassault Systemes, 2011.

[24] Wang Y, Abe S, Latham S, Mora P. Implementation of particle-scale rotation in the 3D-lattice solid model. Pure Appl Geophys 2006;163:1769-85.

[25] Whitney JM, Nuismer RJ. Stress fracture criteria for laminated composites containing stress concentrations. J Compos Mater 1974;8:253-65.

[26] Ippolito M, Mattoni A, Pugno N, Colombo L Failure strength of brittle materials containing nanovoids. Phys Rev B 2007;75:224110-1-0-7.

[27] Landis EN, Nagy EN, Keane DT. Microstructure and fracture in three dimensions. Engng Fract Mech 2003;70:911-25.

[28] Kobayashi M, Tang S, Miura S, Iwabuchi K, Oomori S, Fujiki H. Ultrasonic nondestructive material evaluation method and study of texture and cross slip effects under simple and pure shear states. Int J Plast 2003;19:771-804.

[29] De Proft K, Wells GN, Sluys LJ, De Wilde WP. Combined experimental-computational study to discrete fracture of brittle materials. HERON 2002;47:223-41.

[30] Grosse CU, Finck F. Quantitative evaluation of fracture processes in concrete using signal-based acoustic emission techniques. Cem Concr Compos 2006;28:330-6.

[31] Vidya Sagar R, Raghu Prasad BK, Karihaloo BL. Verification of the applicability of lattice model to concrete fracture by AE study. Int J Fract 2010;161:121-9.

[32] Van Mier JGM. Multiaxial strain-softening of concrete. Part I: Fracture. Mater Struct 1986;19:179-90.

[33] Torrenti jM, Benaija EH, Boulay C. Influence of boundary conditions on strain softening in concrete compression test. J Engng Mech 1993;119:2369-84.

[34] Jansen DC, Shah SP. Effect of length on compressive strain softening of concrete. J Engng Mech 1997;123:25-35.

[35] Van Mier JGM. Multiaxial strain-softening of concrete. Part II: Load histories. Mater Struct 1986;19:190-9.

[36] Van Geel HJGM. Concrete behaviour in multiaxial compression: experimental research. Eindhoven: Technische Universiteit Eindhoven; 1998. p. 1169.

[37] Van Mier JGM, Man H-K. Some notes on microcracking, softening, localization, and size fffects. Int J Damage Mech 2009;18:283-309.

[38] Pramono E, Willam K. Fracture energy-based plasticity formulation of plain concrete. J Engng Mech 1989;115:1183-203.

[39] Etse G, Willam K. A fracture energy formulation for inelastic behavior of plain concrete. J Engng Mech 1994;120:1983-2011.

[40] Menetrey P, Willam KJ. Triaxial failure criterion for concrete and its generalization. ACI Struct J 1995;92:311-8.

[41] Kang HD, Willam KJ. Localization characteristics of triaxial concrete model. J Engng Mech 1999;125:941-50.

[42] Grassl P, Lundgren K, Gylltoft K. Concrete in compression: a plasticity theory with a novel hardening law. Int J Solids Struct 2002;39:5205-23.

[43] Li T, Crouch R. A C-2 plasticity model for structural concrete. Comput Struct 2010;88:1322-32.

[44] Pijaudier-Cabot G, Bazant ZP. Nonlocal damage theory. J Engng Mech 1987;113:1512-33.

[45] Mazars J, Pijaudier-Cabot G. Continuum damage theory - application to concrete. J Engng Mech 1989;115:345-65.

[46] Ju JW. On energy-based coupled elastoplastic damage theories: constitutive modelling and computational aspects. Int J Solids Struct 1989;25:803-33.

[47] Grassl P, Jirasek M. Damage-plastic model for concrete failure. Int J Solids Struct 2006;43:7166-96.

[48] Cervenka J, Papanikolaou V. Three dimensional combined fracture - plastic material model for concrete. Int J Plast 2008;24:2192-220.