Lattice model evaluation of progressive failure in disordered particle composites

Engineering Fracture Mechanics Vol. 57, No. 2/3, pp. 301-318, 1997

~

Pergamon PII: S0013-7944(97)00011-8

LATTICE MODEL FAILURE

EVALUATION

IN DISORDERED

© 1997 ElsevierScienceLtd. All rights reserved Printed in Great Britain oo13-7944/97 $17.oo+ o.oo

OF PROGRESSIVE

PARTICLE

COMPOSITES

B. CHIAIA,t A. VERVUURT and J. G. M. VAN MIER Delft University of Technology, Stevin Laboratory, 2600 GA Delft, The Netherlands Abstract--Results from experimental and numerical investigation on the fractal properties of particle composites are presented. Splitting tests on concrete were carried out and the fractal dimension of the crack patterns detected on microscope images was measured. Numerical simulations of the tests were performed by means of a lattice model, and non-integer dimensions were measured on the predicted lattice damage patterns. The ability of the model to reproduce realistic statistical interactions and self-organization in the propagation of the cracks is discussed. In particular, an increasing Hausdorff dimension was measured during damage development in the lattice network. It is concluded that, next to the stable crack growth provided by fractality in disordered materials, an optimal choice of the percentages of weak and strong microstructural elements, together with their particle-like distribution, may lead to improved mechanical performance of the considered materials. © 1997 Elsevier Science Ltd

1. INTRODUCTION:DISORDER AND FRACTALS THEREEXISTSa class of materials that are heterogeneous at the meso-scale, exhibit strain-softening behaviour and fail progressively by damage-localization and breakdown of their components. Ceramic composites, fibrous materials, concrete, mortar and rocks can be considered to belong to this class of materials. Their global mechanical behaviour is remarkably similar, despite the fact that the mechanisms of failure at the microstructural level are different. Moreover, they are all characterized by a large scatter of experimental results, by considerable size-effects and by a marked difficulty in the reproducibility of the tests. The application of fractal geometry to cementitious materials is relatively new in comparison with the application to rocks and metals. Nevertheless, it has proven very appropriate, due to the multi-scale heterogeneity of these materials [1]. This is reflected in the hierarchical process of failure, in the sense that cracks usually originate from pre-existing micro-defects in the matrix or at the matrix/aggregate bond (micro-level), then develop at the scale of the matrix, involving debonding or even breaking of the closest aggregates (meso-level), and finally extend to the structural scale (macro-level), resulting in the typically disordered (zig-zag) pattern with multiple cracking randomly diffused at the interfaces, through the matrix and through the aggregates. Peculiar toughening mechanisms come into play, like bridging, aggregate interlock and creep [2]. A strong similarity exists between the turbulence in a fluid flow and the fracture in a solid material. The characteristic features of critical phenomena, such as hierarchy of interacting defects, stochasticity and far-from equilibrium instabilities, are present in both processes. The irregularity of the crack trajectories is a well documented and discussed phenomenon [1]. The two potential sources of disordered damage patterns are the intrinsic instability of the fracture process and the effect of the (random) pre-existing material inhomogeneities. An extensive documentation exists on the stability of crack trajectories to small perturbations and it is nowadays widely demonstrated that the quenched disorder of brittle fracture patterns cannot be ascribed to crack-tip instabilities. It can be concluded that the main source of randomness and chaotic behaviour is the inherent disorder of the material microstructure. The overall mechanical behaviour of heterogeneous materials is highly affected by the interplay of disorder and correlations. Any proposed model should not neglect these fundamental aspects. In this paper we address the above issues in more detail, building upon experience with lattice modelling and applications for several years. The paper is not a rigorous and detailed account of the fundamental aspects of lattice modelling, but rather is an attempt to better understand the tAlso at Politecnico di Torino, Department of Structural Engineering, 10124, Torino, Italy. 301

302

B. CHIAIA et al.

mechanical (fracture) behaviour of complex disordered materials in a structure. Many issues and artifacts of lattice modelling have not been completely resolved to date. Nevertheless, the combined experimental/numerical approach advocated in this paper shows clear directions for materials engineering in the future. 2. LATTICE TYPE FRACTURE MODELS

2.1. Outline of the model: mesh generation and implementation of disorder Lattice models have been used over the last 40 years, for various purposes and in many different versions. Hrennikoff [3] introduced a regular triangular lattice of truss elements to solve classical problems of elasticity. However, due to the lack of sufficient computational power, the model remained theoretical. In the late 80s, lattice models were re-introduced by theoretical physicists [4], to calculate problems of conductivity and to model brittle failure in disordered materials. Schlangen and Van Mier [5] were the first to apply lattice models for simulating progressive failure in concrete. Beam elements were used instead of trusses. Due to the fixed connection at the nodes, the equations of a 2D-network of beam elements represent a straightforward discretization of a 2D Cosserat elastic medium, since micro-couples are included. It has been shown by Schlangen [6] that local rotations at the nodes play a fundamental role in the simulated fracture process. Therefore the beam discretization prevails over a network of truss elements. On the other hand, an artificial internal length is imposed by the model, namely the beam length 1, as will be explained in Section 3. When a heterogeneous material is simulated, disorder must be implemented. Several methods can be distinguished. When a regular triangular lattice of beams [Fig. l(a)] is used, heterogeneity must be obtained by means of random distributions of the properties of the beams, since a regular lattice acts homogeneously by nature. However, next to varying the properties of the beams, also a lattice can be generated where an additional parameter of disorder is included, for example a variable beam length [Fig. l(b) and (c)]. The generation of the random lattice is illustrated in Fig. l(c). First, a regular square grid with cell-size s is chosen. In each cell a node is selected at random. Subsequently, the random lattice is generated by connecting the three nodes which are the closest to each other. The connectivity of the nodes is determined by the Voronoi construction [7]. Because of the regularity of the square grid, only a limited amount of neighbouring nodes has to be considered in the numerical procedure. Furthermore, the degree of randomness A [as indicated in Fig. l(c)] can be varied such that, in the limit case where A = zero, a regular square lattice with crossing diagonals is recovered. This case is however never used in the analyses. For A = s, the selected node can be situated in a complete cell of the regular grid. Varying the randomness, the average beam length in the mesh is affected, as illustrated in Fig. 2. The randomness was varied from 0.001s s

\v/

\v /

B

m

\v/

_

I

IA~'~'-"7A~-'--'--7A~"

(a)

(b)

(c)

Fig. 1. Triangular lattices: (a) regular and (b) random. In part (c) the construction of the random lattice with variable randomness is illustrated.

303

Progressive failure in disordered particle composites ;'..;' :

(a) 1.0

%

(b) A=0.Ss

J

.z 0.8

i

--

(c)

o.6 0.4

0.2 0.0 1.100

A=0.Ss

3 L

\ \\

(d) A=0.2s

1.125

1.150

b e a m l e n g t h I/s

Fig. 2. Influence of the randomness A on the average length of the beam elements[16]. Meshes with different randomness are shown in (b), (c), (d).

to 1.0s with steps of 0.1s, and for each randomness at least 175 simulations were performed. The average length of the beams in a mesh of 50 x 50 nodes was calculated. Each simulation is indicated by a black dot in the curve, whereas the average of all the simulations corresponding to each randomness is indicated by the white dots. The solid curve shows the best fit. For the minimum randomness (A = 0.001s), the average length of the beams can be calculated as: l~vg=

1 + ~--

~ 1.14s.

(1)

When the randomness increases, the average beam length decreases, because only the closest nodes are connected. However, when the randomness exceeds 0.5s [Fig. 2(c)], connectivities are not only realized with the direct neighbours, but also with nodes in cells further away. This leads to an increase of the average beam length to lavg~ 1.14S for A -- s. Moreover, it is clear that the scatter globally increases with increasing randomness, as indicated by the individual results of the simulations [the black dots in Fig. 2(a)]. When all beams have the same length, like in the regular triangular lattice [Fig. l(a)], or when the randomness A in the random lattice is chosen equal to zero, no heterogeneity is included in the mesh, and has to be implemented separately. A rather simple way to achieve this is to assign to the beams random strength and stiffness values [4]. This can be obtained by means of statistical distributions (Gaussian or Weibull distributions will be investigated in this paper). Looking at the structure of concrete, it seems more realistic to use a (computer-generated or image-scanned) grain structure as an overlay on top of the lattice [8]. Subsequently, the strength and stiffness of each beam can be determined from the position of the beam in the particle structure, as indicated in Fig. 3(b) and (c). In this figure, " A " represents the aggregate phase, " M " the matrix and " B " the bond between matrix and aggregate. In the case of a computer-generated grain structure, circles are generated according to a Fuller curve distribution derived by Walraven [9]. The particles are placed in a 2D box from which the required grain structure is selected. Two different methods for placing the particles in the box can be applied. In the first procedure random X- and Y-values are chosen for the center of the particle and a check is performed on earlier placed particles whether no overlaps appear. This procedure, however, is very time consuming because a particle has to be replaced (and checked again) as soon as an overlap appears. Especially when a very dense aggregate structure is generated, this may lead to enormous computational efforts [8]. Therefore, in this paper, a different placing procedure has been applied, which is much more efficient. Out of the generated particles, the particle which is to be dropped is selected at random. Instead of choosing both a random X- and Y-value, only X is chosen at random. The Y-coordinate is determined by dropping the particle from the top of the box [Fig. 3(a)]. As soon as the particle touches one of the earlier dropped particles, the Y-coordinate is fixed. To make the particle structure more compact, this placement procedure can

304

B. C H I A I A et al.

be repeated for each particle several times in order to determine the lowest position. Because the generated grain structure does not include all particles (the smallest have to be omitted for computational reasons), a minimum distance between two aggregates is prescribed. The minimal distance is determined empirically, as soon as the box is filled appropriately. Next to the minimum aggregate size, also the step size in the Fuller discretization influences this minimum distance. Two grain structures were used for the simulations presented in this paper. The first contains particles whose diameter is comprised between 0.25 and 8 mm. The second contains particles between 0.25 and 16 mm. In both cases the step size was 0.25 mm, and each particle was dropped l0 times from random positions at the top of the box. Because the beam length has to be at least two or three times smaller than the minimum aggregate size [Fig. 3(c)], the smallest particles (0.25-2 mm) were excluded. When a grain structure with a larger minimum aggregate size is required, the smaller particles can simply be omitted. 2.2. Parameters related to the global elastic behaviour 2.2. I. Regular lattice. The advantage of the lattice model is that it is a simple and transparent model. The parameters needed in the model can be divided into two groups, i.e. those assigned according to the elastic behaviour of the material (theoretical parameters) and those needed in the fracture law of the model (empirical parameters). To describe the global elastic behaviour of each phase in the material, the Young's modulus (E) and the Poisson's ratio (v) of the material are available as input. They have to correspond to the global behaviour of the mesh and can be adjusted by changing the geometrical properties (in-plane height h and out-of-plane thickness t) and the Young's modulus (Eb~am) of the beams. In the lattices presented in this paper, all beams in the mesh have the same geometrical properties and the thickness of the beams is kept equal to the specimen thickness. For the regular lattice, the beam properties can be determined in a rather straightforward manner, since Poisson's ratio is directly related to the height over length ratio of the beams [6]. The beam length l was set equal to l mm in the case of the uniaxial tension simulations and to 1.67 mm in the case of the splitting simulations. Subsequently, a height over length ratio h/l = 1.13 was used, in order to obtain a realistic Poisson's ratio. Heterogeneity in the regular lattice was obtained by varying the ratios of the strengths and of the elastic moduli of the beams in the three respective phases (matrix, aggregates and interfaces), as explained in the previous section. The strength ratios will be discussed in Section 2.3. In this paper, the stiffness ratios have been chosen as follows: EA/EB = 8 and EA/EM = 2, which are considered representative of the most common concrete mixtures (with dense natural aggregates). Subsequently, the Young's modulus of the global mesh can be changed, by varying the local Young's modulus of the beams since the ratio E/Eb~amis constant. The value of the ratios, however, must be maintained.

q

matrix (M)

e J ~ r x / v v x / v ~ ,

x x

~,

(a)

(b)

w_J

(c)

Fig. 3. Particleoverlay(a), projection on top of a regular triangular lattice (b), and assigningproperties to the beams in each phase (c)[8l.

Progressive failure in disordered particle composites

305

2.2.2. Random lattice. Since the beam length is not fixed in the case of the random lattice (Fig. 2), determining the beam height h and Young's modulus Ebeam is more difficult. Therefore, a statistical approach has been followed in which the cell-size s is used to fit the local beam properties in the random lattice (Fig. 4). From the simulations used for determining the average beam length (Fig. 2), also a linear elastic analysis was performed in which one boundary was subjected to a uniform deformation. For each mesh, the height over cell-size ratio his was changed between 0 and 1 with a step-size of his = 0.2. It is mentioned that also ratios up to his = 20 were tested. Because the Poisson's ratio for his > 1.0 converges to values below 0.1, only a few simulations were performed for that range, since these values are not representative for concrete. As well as the variation of h/s, also the local Young's modulus of the beams (Ebeam)was varied. From the linear elastic analysis, the global Young's modulus E and the Poisson's ratio v for the mesh can be calculated. Because the analysis performed on each square mesh is linear elastic, the ratio E/Eb~,m remains constant as soon as h and s are fixed. The results for the Poisson's ratio are given in Fig. 4(a), whereas Fig. 4(b) comprises the results for the Young's modulus. The curves in Fig. 4(a) represent the best fit curves through the results. From Fig. 4(a) it can be seen that there is a considerable influence of the randomness on the measured Poisson's ratio. For increasing randomness in the nodal structure of the mesh, a higher Poisson's ratio is observed. This behaviour seems logical if it is considered that a regular square mesh has a Poisson's ratio equal to zero [3]. However, for h/s ratios larger than approximately 0.Ss, this trend is opposite. In the case of concrete, a global Poisson's ratio of about 0.2 has to be provided. The consequence is that only lattices with a randomness larger than A = 0.25s can be used. From the his ratios and from the global Young's modulus, the beam elastic modulus can be determined from Fig. 4(b). It can be seen that, in contrast with the case of the Poisson's ratio, no influence of A on the ratio E/Eb~amis observed. Moreover, the scatter in the results is small compared to the results presented in Fig. 2. It can be concluded that the elastic modulus of a material is less sensitive to microstructural disorder, with respect to other mechanical quantities (e.g. tensile strength and fracture energy). This is because of localization, which strongly affects strength and toughness and only slightly influences the elastic modulus. Once the global Young's modulus is fixed, the Young's modulus for the beams is fixed as well. However, if another modulus is required, this can be adjusted afterwards, because of linearity of the simulations. 2.3. Fracture law parameters and simulation of material failure Next to the parameters related to the global elastic behaviour of the mesh, parameters related to the adopted fracture law are required. First, the strength values have to be assigned to the beams, based on a statistical distribution or according to their position in the particle microstructure. In the case of normal mixtures, the following strengths were determined as the most appropriate: f~A~ = 10 MPa, f~M) = 5 MPa, f~B) = 1.25 MPa. Note that the ratios are significant and not the absolute values. When a lightweight concrete has to be simulated, the strength values are close to one another, since properties of the constituents are almost similar in the mixture. After generating the mesh and assigning the elastic properties to the beams, the load is applied and a linear elastic analysis is performed. Failure is obtained by removing one beam from the mesh each time. At each load step, the element with the highest stress relative to its strength [(o'/ft)max] is eliminated and (a)

(b)

0.4

O
~ d~ L

0.3 ! 0.2

"t~ ~ ~ ~

-A=I.600s " " A--O.TSOs &--"IA--'0.500S IF~-V A--0.250s : : A--0.0Ols

~

i

"

o.1 0.0 0.00

.~ 0.25

5

0.50

h/s

0,75

1.00

,

J,

o 5

o ~

0.00

/ J

_~jw 0.25

i 0.50

0.75

1.00

h/s

Fig. 4. The beam height as a function of (a) the global Poisson's ratio and randomness, and (b) of the Young's moduli of beams. EFM 57,'2 3 ~G

306

B. CHIAIA

et al.

calculation is re-started after updating the stiffness matrix. The stress a in each beam is calculated according to a very simple law:

F ~ = fl

-A + ~

[Mi,Mj[max) W

'

(2)

where F represents the normal force in the beam, A is the cross-sectional area, Mi and M, are the moments acting upon the two nodes of the beam and W = bh2/6 is the sectional modulus of the beam. The non-dimensional parameters to be determined for this fracture law are the factors ~ and 13. The parameter fl can be determined from the globally measured peak stress in a standard (tensile) test, whereas ~ indicates the role of bending in the fracture law. The factor ~ is less unequivocal since it only affects the tail of the softening diagram [5]. For the simulations presented in this paper, is kept constantly equal to 0.005. The parameter fl depends on the size distribution of the aggregate overlay and on the (average) beam length. The sensitivity of a material to flaws or, inversely, its fracture toughness can be adequately investigated by means of these types of lattice models. The stability of crack propagation is directly tested whenever an element breaks. The released potential energy becomes available to create new failures which again absorb energy. If the released energy is sufficient, another beam fails, and so on. The critical point (brittle failure) is reached when a cascade process is generated. In the adopted approach, the stability is investigated in an indirect manner because, after removal of a beam, the calculation is re-started from zero-load. Moreover, a strength criterion [eq. (2)] is adopted for failure of the elements (toughness is considered to be a macroscopic property derived from random local brittleness, whereas the strength f is defined at the beam level). The global stiffness of the specimen decreases monotonically as the elements progressively fail, and therefore its compliance increases at an increasing rate as members of the network fail (Fig. 5). The critical slope in the compliance curve, corresponding to the onset of catastrophic collapse where the cascade process starts, may be assumed as a material toughness parameter since it is less dependent on the specimen geometry with respect to other toughness parameters. Local instabilities (sudden bursts of released energy) are evidenced in the numerical simulations, in the form of steep descending branches in the load-deformation curve (Fig. 5). It can be argued that these load drops correspond to local snap-backs occurring at the material level, which are difficult to detect experimentally. The main reason is that the real fracture process is a 3D phenomenon, whereas most simulations are carried out in 2D. The higher the strength of the aggregates (as, for example, in the case of phosphorous slag aggregates), the more pronounced are the instabilities. 3. STATISTICAL ASPECTS OF THE LATTICE M O D E L An effective meso-level description of the damage evolution process is provided by the applied lattice model. The cooperative aspects of the phenomenon are also included. Because the model contains the basic features of statistical physics, it can be used to extract information on correlations, criticality and fractality in the fracture process. How can we know in advance if a system is fractal or multifractal? Whenever entropy wins over energy, the resulting structure will be dominated more by randomness than by strict Euclidean order and we might expect to find fractal patterns with a scaling symmetry analogous to that of the Von Koch curve [10]. Multifractal patterns prevail in nature when the interplay between two length scales arises, whereas mathematical fractals, lacking any internal length, can be characterized by a unique value of the self-similarity dimension. It has to be pointed out that, in the lattice model, the beam length l introduces an artificial correlation length superposed to the true internal length lch related to the statistical distributions of the mechanical properties. The internal length scale represents the spatial range of significant mechanical interactions among nearby points. It can be defined statistically as the correlation length of the material microstructure, which decreases with increasing straining and progressive localization. It can be argued that the beam length does not affect the overall behaviour of the composite, provided that it is sufficiently small compared to the size of the specimen. On the other hand, its influence is strong at the material level, and the local breaking rules depend on l (as in

Progressive failure in disordered particle composites

307

(a) 0.010

"/I

0.008

0.006

0.004

0.002

0.000

t

i

t

i

t

t

t

I

,

0

I

i

t

i

I

i

|

t

i

*

t

25

*

,

t

|

i

•

i

.

t

50

i

75

number of broken beams (b) 250 2oo

,m

critical

point

150

8 100

50 0

I

0

i

i

L

I

0.005

i

i

T

i

i

0.0 !0

i

i

i

.

.

.

0.015

.

.

I

t

0.020

-

.

.

~.025

displacement Fig. 5. Plot of compliance vs number of broken beams (a), and load-displacement diagram (b) from numerical simulation of a splitting test.

a micropolar medium). It would be desirable to have l < lch, but the computational efforts would be unacceptable. The beam length represents a lower cut-off in the scaling range and thus directly affects the local fractal dimension. The other correlation length that comes into play with propagating damage in the network is the interaction length ~ between the cracks which, inversely, increases as damage proceeds, ideally tending to infinite (in an infinitely large specimen) when a macroscopic crack suddenly generates from the coalescence of microcracks and pores [see eq. (3)]. Percolation represents the simplest model of a disordered system. The progressive damage accumulation in the lattice network is a typical example of a geometrical phase transition, which can be described as a percolation transition. In the simplest approximation, each beam can be broken with probability p or sound with probability 1 - p . At low p, the network is globally sound [Fig. 6(a)] and the broken beams form disconnected clusters (usually along the particle interfaces) whose linear size may be considered proportional to the interaction length 4. On the contrary, at large p, one or more paths between opposite edges of the lattice exist and the network of beams is totally disconnected. Therefore, a threshold concentration Pc must exist in between, where for the first time a single connected path through broken beams develops [Fig. 6(c)]. It can be said that fracture percolates and Pc is called the percolation threshold. The first macrocrack is defined as the infinite cluster, because its size diverges when the size of the lattice increases to infinity. This

308

B. CHIAIA et

(a)

al.

(b)

(c)

Fig. 6. Firststages of damage(a), intermediatestage (b), and percolationthreshold(c) in a latticenetwork subjected to splitting. corresponds to an infinite interaction length among the microcracks (critical point). As p approaches Pc, ¢ increases as: ¢~lP-Pcl

v,

(3)

and the mean number S of beams belonging to a finite cluster also diverges as: S ~ IP - Pcl-~.

(4)

The exponents v and 7 describe the scaling behaviour of two typical quantities at the critical point and are called the critical exponents. They are universal and do not depend on the structural details of the lattice (e.g. square, triangular or random) but only on the Euclidean dimension d of the lattice (d = 2 for two-dimensional lattices). Exact calculations have been carried out for percolation transitions in 2D lattices, yielding v = 1.33 and y = 2.39 [10]. It has also been demonstrated that the (non-integer) critical exponents can be related to the fractal dimension of the damage patterns. A robust feature of random systems is that, regardless of the local definition of the correlation length, the same scaling exponents describe their asymptotic behaviour. The percolation threshold Pc is obviously affected by the distribution of mechanical properties (strength and modulus) of the beams. The probability of a beam to break initially depends only on the statistical distribution of mechanical properties in the lattice. Then it increases with straining and with the accumulation of broken beams in the surroundings. The interactions between the beams favour an ordered (energetic) state, lowering the percolation threshold. Instead, the random distribution of the beam properties favours a disordered (entropic) state which raises Pc. The transition from fixed (quenched) to evolutive (annealed) disorder has to be taken into account [4]. The self-similarity of percolation clusters was first noticed by Stanley [10]. Since ~ is finite below Pc and never reaches the infinite value in finite-sized lattices, self-similarity can be detected only at length scales smaller than ¢. The crossover from the fractal behaviour at small length scales to a homogeneous behaviour at large length scales is a crucial point to be highlighted. This transition has been revealed experimentally on fracture surfaces of disordered materials [11] and has been described by a continuous variation of the fractal dimension named as geometrical multifractality. Such microscopic scaling transition provides analogous transitions in the scaling of the macroscopic mechanical properties, namely the size-effects undergone by strength and toughness at the structural level [1]. 4. EXPERIMENTAL DETECTION OF DAMAGE Splitting tests were carried out on four different types of concrete to make comparisons with the numerical simulations. The mixtures contained river gravel aggregate with different maximum

Progressive failure in disordered particle composites

309

1 mm ,--,

(b)

Fig. 7. Gray-level image of damage stored during a splitting test (b) on normal concrete. Mosaic obtained by pasting together seven high-resolution images (a).

N

(C)

| mm

! mm

(a) i.,~

j?'

L

(d)

E2

H

Fig. 8. High-resolution image (a), thresholded damage pattern (b), and application of the Box-Counting method to the thinned cracks (c), (d).

Progressive failure in disordered particle composites

311

size (2 and 16 mm, named, respectively, Con2 and Con 16), phosphorous-slag aggregates (ConPS) and lightweight (Lytag) aggregates. Details of the tests are given in ref. [11]. In order to obtain information on the progressively propagating microcracks, a small region around the notch tip was automatically scanned during the loading process [Fig. 7(b)]. After breaking of the specimen, mosaics were created by linking the single images belonging to a prescribed crack path [Fig. 7(a)]. It has been shown that the mosaics represent a more global view of the damage. Therefore, they can be used to describe the transition from the fractal disordered patterns detected at the highest magnifications to the Euclidean domains which are usually observed when the scale of observation decreases. In order to extract the topological information about the damage patterns, filtering of the images is necessary. Starting from complex gray-level images [Fig. 8(a)], pixels belonging to the cracks have to be selected and isolated, that is, thresholded [Fig. 8(b)]. This does not turn into a straightforward procedure since discrimination between cracks and dark particles is extremely difficult. Enhancing the contrast in the image helps to distinguish more clearly the crack boundaries from their surrounding. Because a unique gray-level threshold could not be adopted for a single image, selective filtering has been applied to different parts of the image, choosing for each zone the most appropriate threshold value. This usually ranged between 80 and 135, on a gray-level scale between 0 (black) and 255 (white). Minor corrections were sometimes necessary to eliminate isolated black pixels. After thresholding, binary thinning of the cracks was applied in two different ways. The first one yielded skeletonized cracks by averaging black pixels through the width [Fig. 8(c) and (d)]. A serious drawback for this procedure is the loss of many details from the crack boundaries, resulting in smoother thinned patterns, not properly representing the complexity of the energy dissipation which takes place at the crack lips. Therefore, a more effective algorithm was developed, which was capable of skeletonizing the cracks by isolating their boundaries [11]. 5. BOX-COUNTING METHOD Among the different definitions of fractal dimension, the Minkowski-Bouligand dimension represents the best definition for numerical implementation. The Minkowski-cover of a set is obtained by covering the considered domain by means of a collection of regular Euclidean figures, with dimensions equal to the space in which the fractal is embedded. The fractal dimension is obtained by computing the logarithmic density of the measure of these coverings, as their linear size decreases [1]. For computational reasons, square or rectangular grids are always used. This is called the Box-Counting method. It represents the most general deterministic method because it can be applied either to invasive (fracture surfaces) or to lacunar (damaged ligaments) fractals, whereas other methods can only be applied to single fracture profiles [1]. In order to study the complex patterns resulting from the experimental tests and from the numerical simulations, a general-purpose version of the method has been developed for 2D domains. As the linear size e of the covering boxes decreases, their number N increases and their area E decreases. The box-dimension A is obtained as:

,,.0\ log(1/ i)

-- 2 - !im

log

(5)

In practical applications, instead of dealing with a limit, linear regression is carried out in the log N vs log E plot or in the log e vs log E plot. The slopes ~0 and ~ are respectively obtained. From these values, the box-dimension is yielded as A = - ~0 or as A = 2 - 9. Further details of the algorithm are described in ref. [11]. Two stages of the method applied to the thinned experimental cracks are shown in Fig. 8(c) and (d). Application to the damage patterns obtained from the lattice simulations is shown in Fig. 9(b) and (c). Rigorously speaking, the isotropic shrinking of the covering grids implies self-similar scaling of the covering shape. This could be in contrast with the eventually anisotropic (self-affine or multifractal) properties of the fractal domain. Indeed, direction-dependent scaling arises as the magnification factor decreases. That is, looking at the experimental cracks from farther away, self-affine aspects increase. Therefore, direction-dependent rescaling of the grid should be performed. An iterative procedure should be adopted, since the Hurst exponent H of self-affinity

312

B. CHIAIA et al.

is not known a priori [12]. This relatively awkward procedure can be avoided if anisotropic grids are used, for example by means of rectangular boxes, to enhance the statistical fluctuations in one direction. This approach proved particularly effective in the case of the analysis of low-resolution images, where a strong orthotropy was provided between the crack advancement direction and the orthogonal direction. In this way, it was possible to detect the local fractal dimension in the limit of the smallest grids, as well as the multifractal transition to homogeneous behaviour. In the case of the high-resolution images, because the crack advancement direction is significantly inclined and curved, the orthotropy directions continuously change. Further, the adopted resolution allows for the detection of crack branching and secondary cracks [Fig. 8(a)]. Therefore, it seemed more convenient to use isotropic square coverings of the high-resolution cracks, since damage at that scale shows self-similar rather than self-affine scaling. Analogously, a distinction has to be made for the case of the numerical damage patterns. Owing to the relatively low resolution of the lattice (beam length = 1.00-1.66 mm), a relevant cut-off is present in the scaling regime and fractality can be detected only in a narrow range. Self-affinity is provided by single advancing macrocracks which yield very low fractal dimensions. Self-affine global patterns arise also in the presence of almost homogeneous microstructures where highly localized fracture develops [Fig. 10(a)]. On the contrary, as the degree of disorder in the microstructure increases, diffused microcracking occurs especially at the particle/matrix interfaces and self-similar damage spreads, resembling two-dimensional random walks [Fig. 10(b)]. Percolation-like fracture develops, which provides, at the structural level, increase of both strength and toughness.

(a)

iiiii~ii!iii!;iiill

(b)

Lr

el

J

(c)

Fig. 9. Application of the Box-Counting method (b), (c) to the damage patterns extracted from lattice simulations (a).

Progressivefailure in disordered particle composites

(a)

313

(b)

Fig. lO. Localizedfracture in a homogeneouslattice (a) and percolation-like damage in a disordered lattice (b). Tensile splitting test.

6. FRACTAL ANALYSIS OF THE EXPERIMENTAL CRACK PATTERNS The optical detection of the damage permits to highlight some morphological differences in the fracture behaviour of the four different concretes corresponding to the measured mechanical differences [11]. The interface between large aggregates and matrix represents the weakest link of the four mixtures. In the case of the ConPS, the propagation of the cracks is almost totally controlled by the bond strength and the irregular shape of the aggregates contributes to raise the toughness characteristics. When the aggregate strength is comparable to the strength of the matrix, as in the case of Lytag, interface failure plays a less important role. In this case, a lead role is played by the porosity of the particles, which enhances the bonding strength. Peculiar aspects of a more homogeneous behaviour (at the specimen level) are obtained in the case of Lytag and Con2. Cracks occur through the particles in the Lytag and the fracture patterns in the Con2 appear relatively smooth owing to the small size of the aggregates. On the other hand, microscopic bridging was detected at the scale of sand grains in the Con2, as well as rough fractal cracks were found to develop through the Lytag particles. This means that disorder is equally present in the latter materials, but at a lower level, compared to Con 16 and ConPS. These aspects, unfortunately, cannot be revealed in the numerical simulations, owing to the coarseness of the lattice mesh. The Box-Counting analysis has been performed on the skeletonized images. Non-integer dimensions confirm that, at the proper scales, the phenomenon of fracture in heterogeneous media possesses a well-defined fractal character. A minimum box size equal to one pixel (which means roughly 1.4/~m and 3.7/~m, respectively for the highest and the lowest magnification factor) has been adopted for analysing the thinned cracks. When self-similar damage was detected, through-the-width skeletonizing was often compared with boundary skeletonizing, yielding approximately the same results. Boundary thinning was always preferred in the presence of self-affine patterns, because it allowed to enhance the crack lip's irregular fluctuations orthogonally to the advancement direction. Thereby, the local fractal dimension was calculated. The Box-Counting analysis was also applied to some of the mosaics, whose resolution depended on the number of pasted images. A general trend has been revealed in the four different mixtures: as the resolution of the image decreases, the fractal dimension decreases as well. This shows the transition from fractal microscopic topologies towards homogeneous patterns at larger scales (geometrical multifractality) that occurs in the evolution of damage in the four concretes. Multifractality implies the progressive homogenization of the microstructural fluctuations, which is reflected by the progressive vanishing of the scaling properties of disordered materials for increasing structural size [1].

314

B. CHIAIA et al.

In any case, if the mean values of the measured fractal dimension A are compared (Table 1), it can be concluded that homogenization of the disordered microstructure occurs earlier for the case of light-weight concrete and mortar. Therefore, these mixtures can be characterized by a smaller internal length scale. This may explain the lower scatter which has been found in the measurement of their mechanical properties [11]. 7. FRACTAL ANALYSIS OF THE NUMERICAL CRACK PATTERNS For the analysis presented in this section, a lattice mesh of 3393 nodes was used, with three degrees of freedom per node. The area where no cracks were expected was modelled using conventional four-noded plane stress elements (2 d.o.f, per node) that were available in the finite element package that was used. All analyses were carried out on a Silicon Graphics R4000 workstation (32 Mb main memory). The multifractality of the stress distribution in the beams of a lattice at the critical point has been investigated by theoretical physicists [4]. It is worth clarifying that, in that case, the definition of multifractality was a statistical one, whereas geometrical multifractality represents a topological concept. On the other hand, both approaches have in common an infinite spectrum of (fractional) exponents needed to characterize the scaling of some properties over the whole scale range. Thus, fractal behaviour in the fracture phenomenon is characterized by chaotic damage patterns as well as by the complexity of the stress-strain fields. The relatively coarse mesh (beam size l = 1.00-1.66 ram) introduces a lower cut-off in the scaling behaviour which is about three orders of magnitude larger than the pixel's resolution in the images. Therefore, comparison between experimental and numerical results has to be made with care. In particular, it is not possible to measure the fractal dimension of the final separation surface in the lattice, because of the coarseness of the domain. Also, the fractality of a localized propagating crack is not revealed. Instead, very interesting results come from the fractal analysis of the global damage patterns in the lattice network, which is directly related to the total dissipated energy (Fig. 9). The fractal dimension of the evolutive set obtained collecting the broken beams will be called hereafter simply the damage dimension Ad. This means that the fractal analysis was not used for assessing the geometry of a single macroscopic crack, but rather the complete microcrack pattern, showing the interactive nature of the fracture process in disordered materials, was examined. Globally, the same multifractal trend detected on the experimental images is obtained when looking at the lattice at progressively larger distances (i.e. with decreasing resolution). On the other hand, the beam length represents a coarser cut-off scale with respect to the pixel size, and multifractality is less pronounced in the lattice patterns than in the experimental ones. In the case of the numerical damage patterns, a unique dimension Ad is yielded by the Box-Counting procedure. Thus, self-similar damage comes into play in the lattice model. The fractal characteristics of the damage domains resulting from the lattice computations have been extensively studied. Uniaxial tensile tests (with fixed or rotating boundaries), splitting tests and compression tests have been simulated. Disorder in the microstructure has been implemented in two different ways. In the first approach, normal and Weibull distributions, with and without truncations, and with different variances, have been adopted to distribute the strength and/or elastic moduli of the beams. Because the geometrical heterogeneity is believed to be a fundamental aspect to be modelled, computer-generated particle microstructures were superimposed over the lattice, and mechanical properties were assigned to the beams according to their position in the network

Table 1. Measured fractal dimensions on digitized images Image--,

High resolution (0.71 pix//~m)

Standard resolution (0.278 pix/#m)

Mosaics high resolution

Mosaics standard resolution

Concrete

Min-max

Mean

Min-max

Mean

Min-max

Mean

Min-max

Mean

Con 16 Con2 Lytag ConPS

1.16-1.33 1.08-1.30 1.13-1.35 1.14-1.29

1.25 i.18 1.24 1.18

1.09-1.22 1.07-1.19 1.05-1.20 1.08-1.20

1.14 1.10 1.10 1.14

1.05-1.20 1.05-1.16 1.04-1.15 1.07-1.18

1.13 1.10 1.09 1.13

1.06-1.13 1.02-1.04 1.03-1.08 1.05-1.14

1.08 1.03 1.04 1.07

Progressive failurein disorderedparticle composites 3.0

~"~ r,,["~?(":"~ .<-....~ •

315

......... y = 2.0788 + -1.0027x . . . . Y = 2.4277 + - 1.2593x

2.5 [ |

.~

"\ ~~ ....~ . . . / ~ . . , . ~ ~ ~ ~ .

- + y = 2.6848 + -1.4470X . . . . . y = 2.7419 + -1.4904x

0

+-'+'++

,o

0.5

,

-0.2

,

,

,

,

,

,

0.0

n

.

0.2

.

.

.

.

.

0.4

.

.

.

.

.

.

0.6

.

.

.

.

.

.

.

o.g

.

.

l.O

.

.

1.2

1.4

log(box size) Fig. 11. Log-log plots of number of filledboxes vs linear size of the boxes at various stages of damage. The absolute value of the slope is the fractal dimensionAd.

(see Section 2.1). Beams belonging to aggregates were characterized by the higher stiffness and strength, whereas the lowest properties were assigned to the beams positioned at the particle/matrix interface. The matrix was characterized by intermediate values of strength and stiffness (see Section 2.3). An interesting trend has been revealed, namely the damage dimension Ad increases with the number of broken beams, that is, as damage propagates. In the log N vs log e diagram one observes a progressively increasing intercept as well as a progressively increasing slope of the best-fitting line (Fig. 11). Whilst the first parameter is dependent on the "mass" of the fractal set (i.e. on the number of broken beams), the absolute value of the slope is the fractal dimension of the set, or the damage dimension Ad. Therefore, it can be concluded that the damage locus becomes progressively more disordered with the development of (micro-)cracks. On the other hand, the rate of growth of fractality is maximal at the beginning of the loading process (when entropy prevails over energy) and progressively decreases afterwards, tending to a horizontal asymptote in the post-peak regime (Fig. 12). This matches the typical features of self-organized processes [13]. Chaotic damage accumulation at the first stages leaves place for the self-organization of 2.0

/

,~ 1.5

Bmwnian disorder

1.0 ./: "ill8" "

..... tmi~tiJlUai~¢ti°amlm~itm (fix~) tensi°n . -uni~t.ialtemioa(rol~inll) ...... s p l l t t l n $ ~

-

0.5 otic

~ 0.0

0

self

percolation 500 1000 Number of broken beams

1500

Fig. 12. Evolution of the damage dimension Ad vs the number of failed beams N in different test simulations.

316

B. CHIAIAet al.

Table 2. Mean damagedimensionsAo resultingfrom lattice simulations.Values refer to the last stage before complete fracture Test Gaussian WeibuU Random particle Randomparticle geometry distribution of distributionof distribution distribution beam strength beam strength (dmax = 8 mm) (dm.x= 16 mm) Uniaxial tension 1.47 1.45 1.62 1.55 Fixed boundary Uniaxial tension 1.42 1.34 1.55 1.54 Rotating boundary Tensile 1.40 1.35 1.46 1.48 Splitting test Uniaxial 1.55 Compression microcracks at the critical point (self-organized criticality). Afterwards, the fractal dimension increases very slowly (highly correlated stage). Similar results were obtained by Carpinteri and Yang [14] by means of an analytical K~c-approach in a homogeneous material with randomly distributed microcracks. The asymptotic values of the damage dimension of the lattice patterns are shown in Table 2. It is interesting to note that the asymptotic values of the damage dimension Ad are close to 1.5 for all the numerical simulations. This may confirm the hypothesis of Brownian disorder as the upper limit of complexity in dissipative processes [1]. Extrapolating the results to 3D lattices, energy is considered to be dissipated in a 2.5-dimensional domain. This is intermediate between a process where energy is dissipated on a surface only (LEFM) and a fracture process where energy is dissipated in the complete volume of the material under consideration (which is the domain of plasticity and damage mechanics). The depicted trend has been obtained in all the numerical simulations. The strong similarity existing between experimental and numerical damage patterns, river networks and lightning suggests that Brownian patterns arise in nature as a joint consequence of optimality (in the energy dissipation) and randomness (in the field variables). Besides the common trends in the self-organization of damage, some aspects have to be distinguished between the different simulations carried out by means of the lattice model (see Table 2). In the case of uniaxial tensile tests, higher fractal dimensions Ad are measured when fixed boundary conditions are adopted instead of rotating ones. Further, a larger number of beams fail before final separation. This can be explained with a larger correlation among the microcracks in the case of rotating boundaries, where a single main crack develops from one side of the specimen to the other [Fig, 13(b)]. Fixed boundaries allow for the formation of (at least) two macro-cracks, thereby implying a larger degree of disorder in the damage process [Fig. 13(a)]. Earlier experimental tests and numerical simulations by Van Mier e t al. [15] also provided-higher fracture energy values in the case of the fixed boundaries specimens. Splitting simulations yielded, on average, lower damage dimensions Ad. This might be explained by the higher localization of damage which occurs in this geometry (Fig. 10). Obviously, a lower damage dimension Ad does not imply a lower local dimension of the (self-affine) crack lips, which should be independent of the test geometry. Unfortunately, the lower cut-off in the scaling range represented by the beam length does not allow the measurement of the local fractal dimension of the crack lips (where self-affinity comes into play). The maximum aggregate size scarcely influences the values of the damage dimension. Instead, it affects the scale range where fractal patterns arise. It may be concluded that 2.5 represents a universal dimensionality for the energy dissipation domain, regardless of the microstructural details. Increasing the beam size from 1 mm to 1.66 mm does not affect the values of the fractal dimension, but causes a shift of the scale range where fractality emerges. Therefore, if the size of the specimen is enlarged proportionally to the beam size, similar damage patterns are obtained. Nevertheless, if only the beam size is augmented, a larger internal length is artificially imposed to the material and the topologies of damage differ substantially. This means that the effects of disorder on the mechanical properties become more relevant as the internal length increases. If standard statistical distributions are used to assign properties to the beams, the damage dimension Ad increases with their variance up to a certain limit, depending on the total number

Progressive failure in disordered particle composites

(a)

317

(b)

Fig. 13. Different damage patterns in the presence of fixed (a) and rotating (b) boundaries in a uniaxial tensile test.

of beams. Afterwards, Ad remains almost constant. Gaussian distributions yield higher dimensions with respect to Weibull distributions, because too many weak beams imply a more rapid transition to percolation threshold in the network. The damage dimension increases if, instead of classical statistical distributions, a particle microstructure is projected on the lattice and properties are assigned to the beams according to their position on the projection. Thus, the spatial organization of a particulate composite with relatively weak interfaces and tough aggregates seems to be an effective remedy against the intrinsic brittleness of the constituent materials. Indeed, a more homogeneous composite like light-weight concrete should be modelled with an interface strength comparable to the aggregate strength (through-the-particle cracking frequently occurs). In this case, the crack-arresting capacity has to be found into the micro-roughness of the Lytag particles, whose fractal character cannot be revealed from the lattice model due to the coarseness of the mesh. 8. CONCLUSIONS Some interesting suggestions for specific "tailoring" of new high-performance cement-based materials can be drawn from the lattice simulations. First of all, the smoothing of the stress singularity induced by the fractality of propagating cracks provides the stable-crack growth behaviour typical of heterogeneous materials. Moreover, the transition to Euclidean behaviour, and the consequent vanishing of the disorder-induced effects (size-scale effects) is governed by the interplay between internal length and structural size. Raising the strength of all beams in the lattice increases the ultimate load of the specimen but does not decrease its sensitivity to flaws. Instead, selectively increasing the strength of some beams can raise both strength and toughness in the specimen. In particular, the spatial organization in a particulate microstructure seems to play a fundamental role in determining the material characteristics. Increasing the degree of disorder (which could be set equal to the variance of the statistical distributions of the mechanical properties or to the randomness of the particles geometrical distribution) decreases the fracture sensitivity, i.e. a larger number of beams randomly fail prior to the critical development of an unstable crack. The amount of weak beams (e.g. those belonging to the interfaces) has to be set within an optimal range. It must not be too small, in order to ensure sufficient random distributed microcracking, nor too large, because otherwise the weakest chain immediately percolates through the lattice. It can be said that the weakest link approach (Weibull) appears suitable to model the

318

B. CHIAIA et al.

material failure only if a very large percentage of weak b o n d s is present in the material, or if homogeneity prevails. It fails to predict failure in the case of real materials, whose microstructures are heterogeneous, a n d possess several t o u g h e n i n g mechanisms. The results from the fractal analyses are in complete agreement with earlier observations regarding the increase of fracture energy u n d e r fixed loading platens in a uniaxial tensile test, as c o m p a r e d to the case where the two loading platens are free to rotate. This m e a n s that the fractal analysis, where the geometry o f all interacting microcracks in a plane is considered, a n d referred to as d a m a g e d i m e n s i o n in the paper, is a valuable p a r a m e t e r for assessing the t o u g h e n i n g m e c h a n i s m s in disordered materials.

REFERENCES 1. Carpinteri, A. and Chiaia, B., Multifractal nature of concrete fracture surfaces and size effects on nominal fracture energy. Mater. Structures 1995, 28, 435~,43. 2. Van Mier, J. G. M., Mode I fracture of concrete: discontinuous crack growth and crack interface grain bridging. Cement and Concrete Research 1991, 21, 1-15. 3. Hrennikoff, A., Solution of problems of elasticity by the framework method. J. appl. Mech. 1941, 12, 169-175. 4. Herrmann, H. J. and Roux, S., Statistical Models for the Fracture of Disordered Media. Elsevier Science, 1992. 5. Schlangen, E. and Van Mier, J. G. M., Experimental and numerical analysis of micromechanisms of fracture of cement-based composites. Cement and Concrete Composites 1992, 14(2), 105-118. 6. Schlangen, E., Computational aspects of fracture simulations with lattice models. In Fracture Mechanics of Concrete Structures, ed. F. H. Wittmann. AEDIFICATIO, Freiburg, 1995, pp. 913-928. 7. Moukarzel, C. and Herrmann, H. J., A vectorizable random lattice. HLRZ, KFA Jiilich, preprint, 1992. 8. Vervuurt, A., Van Vliet, M. R. A., Van Mier, J. G. M. and Schlangen, E., Simulations of tensile fracture in concrete. In Fracture Mechanics of Concrete Structures, ed. F. H. Wittmann. AEDIFICATIO, Freiburg, 1995, pp. 353-362. 9. Walraven, J. C., Aggregate interlock: a theoretical and experimental analysis. Ph.D. Thesis, Delft University of Technology, The Netherlands, 1980. 10. Stanley, H. E. and Ostrowsky, N., On Growth and Form: Fractal and Non-Fractal Patterns in Physics. Martinus Nijhoff, Dordrecht, Germany, 1985. 11. Vervuurt, A., Chiaia, B. and Van Mier, J. G. M., Damage evolution in different types of concrete by means of splitting tests HERON, 1995, 40(4), 285-311. 12. Mandelbrot, B. B., The Fractal Geometry of Nature. W. H. Freemann and Company, New York, U.S.A., 1982. 13. Bak, P., Tang, C. and Wiesenfeld, K., Self-organizedcriticality: an explanation of 1If noise. Phys. Rev. Lett. 1987, 59, 381-384. 14. Carpinteri, A. and Yang, G. P., Fractal dimension evolution of microcrack net in disordered materials. Theor. appl. Fracture Mechanics, to be published, 1996. 15. Van Mier, J. G. M., Vervuurt, A. and Schlangen, E., Boundary and size effects in uniaxial tensile tests: a numerical and experimental study. Fracture and Damage of Quasibrittle Structures, ed. Z. P. Ba~ant, Z. Bittnar, M. Jirb.sek and J. Mazars. E&FN Spon, London/New York, 1994, pp. 289-302. 16. Vervuurt, A., Interface fracture in concrete. Ph.D. Thesis, Delft University of Technology, 1997 (Received 6 February 1996)