- Email: [email protected]
3D lattice type fracture model for concrete
Engineering Fracture Mechanics 70 (2003) 927–941 www.elsevier.com/locate/engfracmech
3D lattice type fracture model for concrete G. Lilliu
a,*
, J.G.M. van Mier
a,b
a b
Microlab, Delft University of Technology, Delft, The Netherlands Institute for Building Materials, ETH, CH 8093 Z€urich, Switzerland
Received 21 February 2002; received in revised form 4 July 2002; accepted 17 July 2002
Abstract A 3D beam lattice model is used for simulating fracture processes in concrete, which is schematized as a three-phase material (aggregate, interfacial transition zone and matrix). Numerical experiments are conducted varying the particle density. The obtained results show that, with increasing particle density, the peak-load decreases and the response is more ductile. This appears to be related to the amount of de-bonding. As a matter of fact, when the strength of the interface transition zone is set equal to the strength of the matrix, neither the peak-load nor the ductility of the lattice response are influenced by the particle density. 2002 Elsevier Science Ltd. All rights reserved. Keywords: 3D lattice; Particle density; Peak-load; Ductility; Interfacial transition zone
1. Introduction Fracture processes due to mechanical loading are a major cause of damage in concrete structures. An adequate comprehension of these processes is important not only for ensuring safety of the construction, but also for optimizing the material behaviour, which would apport beneficial economical effects. Apparently, crack formation and propagation in concrete are affected by its highly disordered structure. Therefore, this structure should be taken into account modelling the material behaviour. At present fracture is modelled at different levels of observation, ranging from the macroscopic level, where material is considered homogeneous and non-linear constitutive relations are used, to meso, micro, or even to atomic or sub-atomic scale. At the meso-level of observation concrete is regarded as a two or three-phase material. A model which allows a straightforward implementation of the material heterogeneity at the meso-level is the lattice model. Such type of material modelling has already been used in the 1940s for simulating problems of elasticity [1]. More recently it has been introduced in the field of statistical physics for studying fracture of disordered materials [2]. Lattice models assume usually a linear-elastic material constitutive relation and, for this reason, they are conceptually very simple. Since the lattice models are a simple and straightforward tool for understanding the physics of fracture processes, they are considered particularly appealing to the authors.
*
Corresponding author. E-mail addresses: [email protected] (G. Lilliu), [email protected] (J.G.M. van Mier).
0013-7944/03/$ - see front matter 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 3 - 7 9 4 4 ( 0 2 ) 0 0 1 5 8 - 3
928
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
Such models, at least in their 2D version, present a characteristic brittleness of the response in terms of load–displacement diagrams. Apparently, this is a consequence of the assumed brittleness of the single element of the network which describes the material. However, it is also important to consider carefully to what extent the material is reproduced realistically in the model. In fact, in order to reduce the computational time, the adopted meshes are usually too coarse. The consequence of this is that only the largest fractions of the particle distributions are reproduced (with diameter da P 2 mm), while in real concrete grains with size of the order of micro-meters are present (e.g. cement grains). Furthermore, the adopted particle distributions, usually computer-generated, hardly are as dense as in real concrete. Last but not least, real fracture processes are 3D and not 2D. Presumably, both the presence of small particles and the 3D effect can increase the ductility of the lattice response. In fact, since more branching and bridging in the crack pattern is expected, the amount of energy dissipated in the post-peak region of the load–displacement diagram should increase. In the past many different types of meso-level models have been developed, for example [3–7]. In all these models the particle structure of concrete was included in a finite element mesh. In some of these models the interface between aggregate and matrix was modelled with zero-length spring elements. These elements had a normal softening component and a shear component along the surface of the aggregate [5,6]. The main objection against these models is that they are unsuitable to explain softening, since softening is already included at the meso-level. An alternative approach was that proposed by [4]. Based on linear-elastic fracture mechanics, crack propagation along the interface and the matrix could be described. The disadvantage of this model is that initial cracks must be defined at the beginning of the analysis, since no crack nucleation criterion is included. Lattice type models have been applied by others, too [3,7]. Schorn [3] used truss elements, which require some numerical measures to avoid instability during fracture propagation. Again, softening was introduced at element level in the model of Bolander [7], who applied a beam lattice with a step-wise softening law for the elements. In this paper results are shown of simulations made with a 3D version of the so-called ÔDelft latticeÕ. This is an extension of a previous 2D version [8], which has been extensively used in the past decade to simulate different experiments [9–12]. In this model beam elements are used for avoiding numerical instabilities. At element level, an elastic purely brittle fracture criterion is adopted, which makes the model suitable for studying the effect of material structure on softening. An additional advantage is that the particle structure can be included in a simple and straightforward manner. Since 3D lattice simulations with relatively fine meshes would be practically impossible with the standard software, the model has been implemented in a code for parallel computing [13]. The case considered is a prism subjected to uni-axial tension. The effect of varying particle density is studied for two cases: in the first case, the interfacial transition zone (ITZ) has a lower strength than the matrix; in the second case, the same mechanical properties as the matrix are assigned to the interface. Thus, in the latter case, the material is considered as two-phase material rather than containing three distinct phases.
2. The model The material is schematized with a network of Bernoulli beams. Lattices can be regular or random. In the first case, all elements of the mesh have the same length. In the second case, after the nodes have been positioned randomly inside the cells of a grid with regular cell size s, they are connected by means of the Voronoi construction [14] (see Fig. 1). For each node, a number of neighboring nodes is checked for possible connectivities (36 in the 2D case, 178 in the 3D case). In random lattices, where all elements of the mesh have different lengths and, thus, different stiffnesses, disorder is already introduced at the ‘‘geometrical’’ level. The amount of such disorder can be varied by placing the nodes inside sub-cells of the original cells. Randomness is defined as the ratio between the size of the sub-cell and the main cell, i.e. A=s. The cross-section of the beams is rectangular and circular in the 2D and 3D case, respectively. The height (in the
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
929
Fig. 1. (a) Voronoi construction for a 3 3 2D mesh with A=s ¼ 0:5; (b) generation of the connectivities for one node of the mesh shown in (c); 3 3 3 3D mesh with A=s ¼ 0:5.
2D case) or diameter (in the 3D case) of the cross-section, h, must be chosen properly in order to obtain the same overall PoissonÕs ratio as the real material [11,15]. Heterogeneities can be introduced in a straightforward manner. The mesh is overlaid on top of a computer-generated particle distribution and different mechanical properties are assigned to the beams falling in each phase. Commonly used ratios for the YoungÕs moduli and tensile strengths are Ea =Em ¼ 70 000=25 000, Em =Eb ¼ 25 000=25 000 and ft;a =ft;m ¼ 10=5, ft;m =ft;b ¼ 5=1:25 (where a, m and b refer to aggregate, matrix and bond, respectively). The resulting density of the aggregates (ratio between number of aggregate beams and total number of beams) is usually much lower than the particle density in real concrete. One of the reasons is that the particles are placed randomly, starting from the fraction with largest size, in a way that a minimum mutual distance is respected. As a consequence, placement of the smallest particle fractions becomes more and more difficult when the density increases. Another reason is that part of the particle volume falls inside the interface (bond zone), which is much thicker in the lattice than in reality. Finally, since the adopted meshes are relatively coarse, the smallest particles are automatically discarded (see Fig. 2(b)). In fact, the length of the beams and the minimum grain size, lbeam and da;min , should follow the minimum requirement lbeam 6 da;min =3. However, this would imply an enormous increase of the number of degrees of freedom and of the computational time. In lattice simulations, fracture processes are reproduced by sequential removal of elements from the mesh and it is assumed that such elements have linear-elastic behaviour until failure. At each step, the element to be removed is that with maximum value of the ratio between the effective stress r , and the tensile strength
Fig. 2. (a) Introduction of heterogeneities by overlaying the lattice mesh on top of a computer-generated particle distribution; (b) distinction between aggregate, matrix and bond phase according to the location of a lattice element in the particle structure.
930
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
Fig. 3. Smoothening of a force–displacement diagram obtained from a lattice simulation.
ft . Despite the assumed brittleness of local failure, global structural softening can be obtained. The effective stress is computed at the ends of the lattice beam, according to: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Mx2 þ My2 F ; r ¼ þ a A W where A and W are the area and the section modulus of the cross-section and F, Mx and My are the internal axial force and bending moments. Thus, it is assumed that only mode I fracture occurs locally. The coefficient a can be used to adjust the contribution of the bending moment in a beam. For 2D lattice analyses it was shown that small values of a give a long tail in the stress-deformation response [10]. However, if the value of a is taken equal to zero, peeling effects might occur in the fracture process in case of a rectangular lattice [16]. A value commonly used in 2D lattice simulations is a ¼ 0:005. It can be shown that in the 3D analyses reported in this paper a variation of a in the range 0.00–0.05 does not affect the crack pattern. Therefore, in this paper the value a ¼ 0 will be adopted, which means essentially that a normal stress criterion is used. For sake of simplicity, torsion is not included in the fracture law. Finally, a remark concerning the load–displacement diagrams resulting from lattice analyses should be made. They present a characteristic zig-zag pattern, which shows that local instabilities occur during crack formation and propagation. Nevertheless, the diagrams obtained from lattice simulations, with their drops
Fig. 4. Case study: (a) computer-generated particle distribution; (b) positioning of the specimen inside the volume Vtot .
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
931
of load, are quite similar to those recorded in the laboratory during experiments [17]. However, when a specific control method is adopted for conducting the experiment, stable crack propagation can be obtained. Still, local instabilities may occur, but they give rise to small indents in the curve only. The zig-zag in the lattice load–displacement diagrams is usually overcome by applying a smoothening technique as shown in Fig. 3. For comparison of various results, smoothening is quite helpful. The inconvenience of this procedure is that information about possible snap-backs is lost and that the amount of dissipated energy is overestimated. In the model the irreversible deformations due to geometrical mismatch of the crack faces during cyclic loading [18] are not taken into account. A way for including permanent deformations is to use a phenomenological approach like in [19].
3. Case study As case study a cube subjected to uni-axial tension is analysed (see Fig. 4(b)). The edge of the cube is 24 mm long. For generating the lattice mesh, the size s ¼ 0:75 mm of the gridpcells ffiffiffi was chosen, with randomness A=s ¼ 0. Therefore, the diagonal beams are the longest, with lbeam ¼ 3s ¼ 1:3 mm. This value of the randomness was selected in order to isolate the effect of the disorder induced by the particles from the geometrical disorder. The nodes of the upper and bottom face of the cube are supported in the X, Y, and Z direction, and an uniform displacement in the Y direction is applied to the nodes of the upper face. The choice of the geometry of the lattice was dictated only by the necessity to limit the computational time. In fact, the number of nodes in the mesh is 35 937, resulting in 431 244 degrees of freedom, which could be handled using 32 processors of a CRAY T3E. Obviously, such a small cube cannot be representative of the overall material behaviour, since boundary effects can be expected. Furthermore, only particles with da P 4 mm can be included, which are relatively big in comparison with the size of the specimen. In this case the material is not Ôstatistically homogeneousÕ and its internal structure (e.g. spatial distribution of particles) might play a major role, causing a large scatter of the outcome from the simulations. Nevertheless, such a numerical experiment is useful for a qualitative understanding of 3D fracture processes. An important parameter to be assigned in lattice analyses is the dimension of the cross-section, which influences the overall PoissonÕs ratio. In random lattices the value of m varies with the randomness of the mesh and when the randomness is A=s ¼ 0 never reaches the value 0.2 [11,15], which is the value usually adopted for concrete. However, it was found that m becomes independent from the randomness when h=s 0:8. For this reason, it was decided to adopt the same value h ¼ 0:577 mm like in simulations of 2D regular lattices with beam length lbeam ¼ 1 mm [20]. It should be mentioned that with the adopted geometry of the beams beam theory is not applicable any longer. This fact is ignored in the analyses to come. As a matter of fact the model is a very rough approximation of reality but this appears throughout the entire approach.
4. Fracture process with sparse particle distribution 4.1. Aggregate structure A particle distribution is considered with maximum grain size da;max ¼ 12 mm and volumetric density Pk ¼ 75%, which is realistic for concrete. A Fuller grading curve is adopted as basis for the computergenerated aggregate structure. The number of particles with diameter di1 6 da 6 di is obtained following qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi di =da;max di1 =da;max Vtot ; Pk
932
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
Fig. 5. (a) Grading curves of the computer-generated particle distribution shown in Fig. 4(a); (b) phase fractions after overlaying the particle structure on top of the 3D lattice mesh with randomly varying position of the specimen inside Vtot .
where Vtot is the volume to be filled (in the example Vtot ¼ 72 72 72 mm3 ). In order to discretize the Fuller curve, the range 0 mm 6 da 6 12 mm of diameters is divided into intervals of 0.5 mm, as shown in Fig. 5(a). As an approximation, all grains falling in one interval of the distribution are considered to have the same diameter, namely, the integer value of the interval. As a consequence, the grains falling in the interval 0 mm 6 da 6 0:5 mm are neglected and the resulting volumetric density decreases to Pk;eff ¼ 60%. Next, the grains are positioned randomly inside the volume Vtot . Starting from the biggest diameter, the particles are placed in a way that a mutual minimum distance, 1:1ðda;i þ da;j Þ=2, between two adjacent particles is guaranteed. Since positioning of particles becomes difficult when the smallest fractions are considered they are usually neglected. Finally, geometrical considerations (see paragraph 2) impose a lower bound to the minimum diameter of the particles (4 mm 6 da ), causing a further decrease of density, until the value Pk;eff ¼ 34%. Once the particle distribution has been generated, the specimen can be placed in the volume Vtot (see Figs. 2(a) and 4 for a 2D and a 3D case, respectively) and heterogeneities can be implemented. After the lattice has been overlaid on top of the particles a new density Pk;latt can be defined as the ratio between the number of aggregate beams and the total number of beams. Fig. 5(b) shows how the fractions of the three phases vary with randomly varying the position of the specimen. As already men tioned in paragraph 2, Pk;latt 6 Pk;eff . In paragraph 4.2 the results obtained from a case where Pk;latt ¼ 29% will be reported. 4.2. Results For the case Pk;latt ¼ 29%, the smoothened force–displacement diagram, until a displacement of 25 lm is shown in Fig. 6(a). A detail of the same diagram, until 5 lm displacement, is shown in the inset. The diagram shows two main branches, namely, the pre-peak and the post-peak branch. In the pre-peak regime, after the first beam element has failed (at the stage marked with a, which corresponds to 20% of the peakload), the diagram shows a kink, which is followed by a nearly linear part (b–c). A strong non-linearity (c– d) can be observed before the peak-load. In the post-peak regime the diagram shows subsequently a steep drop of load (until the stage marked with e) and a long tail. In order to enlighten the development of the fracture process the number of beams belonging to matrix or bond phase, which have failed at each loading stage, is shown in Fig. 6(b). Also aggregate beams fail during the fracture process. However, only 2% of the total number of beams failed at 25 lm displacement are aggregate elements. This value is negligible when compared with the percentages of failed matrix and bond elements, which are 27% and 71%, respectively
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
933
Fig. 6. (a) Force–displacement diagram obtained from a lattice simulation with Pk;latt ¼ 29%; (b) cumulative number of failed beams versus displacement.
(see Fig. 14(a)). Firstly, it can be observed that more interface elements rather than matrix elements fail. Although the first matrix beam fails at point b (60% of the peak-load), practically only de-bonding occurs until point c of the diagram is reached (85% of the peak-load). As a matter of fact, the reason for debonding to prevail is twofold: • stress concentrations occur next to the particles due to the difference in stiffness between the aggregate and the bond/matrix phases, and, • the interface has the lowest strength in the three-phase material. After point c, cracks start to appear in the matrix and propagate rapidly after the peak-load. This corresponds to a steep drop in load, directly following the maximum load. The same phenomenon can be observed in the laboratory and corresponds to macroscopic crack propagation. Finally, prevailing failure of matrix elements is detected along the tail of the softening branch of the diagram. For visualizing the crack process, the crack patterns at stages d, e and at the value 25 lm are shown in Fig. 7. There, an exploded view of the specimen, after it has been cut into 0.75 mm thick slices, is given. The failed beams, falling inside each slice, are projected on a plane, and represented in white or black, depending if they are interface or matrix beams, respectively. From Fig. 7 three stages can be clearly distinguished: • first stage: cracking is due mostly to distributed de-bonding, • second stage: major cracks have formed, • third stage: no relevant changes in the fracture process occur anymore, but crack localization becomes more pronounced.
5. Effect of the particle density on the fracture process 5.1. Aggregate structure As already mentioned in the previous paragraph, high packing densities cannot be obtained with random positioning of the grains (unless much time consuming methods are used). Although yet not as high as in reality, the density of the computer-generated particle distributions can be increased reducing gradually the volume containing the initial sparse packing, by means of dynamic impact models. Using the program
934
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
Fig. 7. Exploded view of the crack patterns in the specimen with Pk;latt ¼ 29% at different loading stages: (a) at peak-load; (b) at 4 lm displacement (point e in the force–displacement diagram); (c) at 25 lm displacement (end of the tail in Fig. 6(a)). The failed bond and matrix elements are represented in white and black, respectively.
SPACE [21], which follows such an approach, it was possible to generate particle distributions with Pk;eff ¼ 48% and 62%, respectively. The difference with the particle distribution mentioned in paragraph 4 is that here the diameter of the particles varies continuously in the range 4 mm 6 da 6 12 mm. Fig. 8(a) shows the number of particles belonging to each fraction of the sieve curve, for the three different densities. In Fig. 8(b) the variation of the phase fractions with the density are shown. For computing the number of beams falling in each phase, one position of the specimen inside the volume containing the particles has been adopted (the same position like that considered in the previous paragraph). Once the lattice is overlaid on
Fig. 8. (a) Number of particles for three different values of the volumetric aggregate density, Pk;eff ; (b) phase fractions of the 3D lattice . as function of Pk;eff
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
935
top of the particle distribution, the densities Pk;latt ¼ 39% and 53%, corresponding respectively to Pk;eff ¼ 48% and 62%, are obtained. The matrix phase decreases with increasing particle density, while the bond phase increases, although slower than the aggregate fraction. The two densities, Pk;latt ¼ 39% and 53%, correspond to the two extreme situations of comparable amount of aggregate and matrix phase, and comparable amount of matrix and bond phase, respectively. 5.2. Results The outcome from the lattice simulations for the values of density Pk;latt ¼ 29%, 39% and 53%, is shown in the diagrams of Fig. 9. The dimensionless load–displacement diagrams, where load and displacement are divided by the values at the peak, are shown in the inset of Fig. 9(a). The number of failed beams versus displacement is shown in Fig. 9(b). Finally, the crack patterns for Pk;latt ¼ 39% and 53%, respectively, are shown in Figs. 10 and 11. From such figures it appears that both the peak-load and the post-peak behaviour are influenced from the particle density. Nevertheless, the fracture process still develops as described in the previous paragraph. Initially, only de-bonding occurs. Next, cracking localizes in the matrix (see Fig. 9(b)) and longer cracks develop in the specimen (see Figs. 7, 10 and 11). Crack localization corresponds to the initial part of the post-peak branch in the load–displacement diagram. A relevant difference among the three cases should be mentioned. For Pk;latt ¼ 29% and 39%, crack localization corresponds to a steep branch after the peak, and matrix failure prevails on interface failure in the tail of the softening curve. For Pk;latt ¼ 53%, after the peak the load decreases gradually while the displacement increases, which produces a more ductile response. Furthermore, interface failure still prevails in the tail of the softening curve. In the next two paragraphs, the effect of the particle density on the peak-load and the ductility of the lattice response, respectively, will be discussed in more detail. 5.2.1. Effect on the peak-load It can be observed from Fig. 9(a) that the peak-load decreases with particle density. This can be explained as follows. The carrying capacity of the specimen depends on the strength of the matrix–bond system. The interface fraction is the weakest component of such a system and contributes to lower its strength. Thus, since the fraction of bond elements in the lattice increases with the density, it can be expected that lower peak-loads are achieved with higher particle densities. Also another factor can contribute
Fig. 9. (a) Force–displacement diagram obtained from lattice simulations with Pk;latt ¼ 29%, 39% and 53%; (b) cumulative number of removed beams versus displacement.
936
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
Fig. 10. Exploded view of the crack patterns in the specimen with Pk;latt ¼ 39% at different loading stages: (a) at peak-load; (b) at 4 lm displacement; (c) at 25 lm displacement.
Fig. 11. Exploded view of the crack patterns in the specimen with Pk;latt ¼ 53% at different loading stages: (a) at peak-load; (b) at 4 lm displacement; (c) at 25 lm displacement.
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
937
to this. At the peak-load the fracture process starts to localize in the matrix. At this stage the structure contains a certain amount of damage, which is related to the amount of de-bonding in the pre-peak regime. As a matter of fact, with increasing number of bond elements more of such elements are removed before cracks start to propagate in the matrix (see Fig. 9(b)). The consequence of this is that the structure is more damaged and its carrying capacity decreases. The dependence of the peak-load on the volume and strength of bond phase will be discussed again in paragraph 6. There it will be shown that, if no bond phase is considered, the peak-load does not depend on the particle density. 5.2.2. Effect on the ductility Here, ductility is defined as the ratio between the energy consumed during the fracture process after and before the peak. The results show that an increase of the particle density produces an up-lift of the softening branch in the dimensionless load–displacement diagram (see inset in Fig. 9) and, thus, an increase of ductility. This is more remarkable for the highest value of the density, Pk;latt ¼ 53%, which approaches the real particle density in concrete (75%). In this case de-bonding plays an important role also in the postpeak region. In this region, not only the number of failing bond elements increases continuously, but even with a higher rate than the matrix elements (see Fig. 9(b)). As pointed out in the introduction, a higher ductility is expected when more branching and bridging occur [22]. Both these phenomena were observed in impregnation experiments during the tail of the softening curve [23,24]. Indeed the current analyses show that more branching, which implies more de-bonding, contribute to increase the ductility of the lattice response. Due to the difficulty in visualizing the crack patterns resulting from the simulations, it is harder to verify that also bridging can be reproduced in lattice simulations. However, matrix failure in the postpeak regime might be related to flexural cracking due to bridging, and not only of macroscopic crack propagation. Thus, realistic particle distributions seem to be essential to obtain a more ductile response. However, some comments should be made. In Fig. 14, the percentage of failed elements falling inside each of the three phases is represented as function of the density (such percentage is computed on the total number of elements removed at a displacement of 25 lm). Fig. 14(a) shows that the percentage of aggregate elements failing during the experiment increases with the density. For Pk;latt ¼ 53%, such percentage reaches the value 4.6%, while the percentages of failed matrix and bond elements are 13.7% and 81.7%, respectively. Although the percentage of failed aggregate beams is negligible when compared with the percentage of failed interface elements, it cannot be ignored when compared with the percentage of failed matrix elements. It is questionable whether the fracture of aggregate beams can affect somehow the ductility of the lattice response. It might be expected that in normal concrete, like is the material considered here, fracture is limited only to interface and matrix. Nevertheless, since also in experiments aggregate failure sometimes occurs, the lattice results still appear quite realistic. The fact that also aggregate elements are fractured might be related to stress concentrations induced by the geometry of the grains. As a matter of fact, since the lattice meshes are always relatively coarse, the contour of the grains is never smooth but shows irregularities (see Fig. 12(b)). Most probably, a more appropriate choice of mechanical properties would have avoided aggregate fracture. For example, it might have been sufficient to assume a higher aggregate strength, i.e. ft;a ¼ 4ft;m rather than ft;a ¼ 2ft;m .
6. Effect of the bond strength on the fracture process From the results obtained with varying particle density it appears that the ITZ, which is the weak layer linking the aggregates to the matrix, plays a major role in the fracture process. The importance of the ITZ is evident also from Fig. 12(b), where a detail of a 2D regular lattice with beam length 0.25 mm is shown. In Fig. 12(b) the particle distribution has a density Pk;eff ¼ 83%, with minimum diameter 1 mm and maximum diameter 12 mm. For such a high density, even though a very fine lattice mesh is used, continuous paths of
938
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
Fig. 12. (a) Schematization of the ITZ in the material with a series model; (b) detail from a lattice mesh where percolation of the bond phase occurs.
bond elements can be detected. In such case the cracks in the bond phase will coalesce, causing formation of macrocracks (percolation). There is quite some debate on the thickness of the real ITZ, which varies in the range between 30 and 50 lm. In the current model the bond zone is a relatively wide zone (at least 10 times as large as the real ITZ), which should be better regarded as a region of average properties among the three phases, namely aggregate, ITZ and matrix. The bond element can be schematized in a simple manner adopting a system of longitudinal springs connected in series (see Fig. 12(a)), with lengths la , lb , and lm , respectively. For such a system, the equivalent elastic modulus Eb;eq , follows l la lb lm ¼ þ þ Eb;eq Ea Eb Em and the equivalent strength is the same as the ITZ strength which is the weakest element of the system. In [11] it was shown that the elastic modulus of the ITZ hardly affects the fracture process, while the strength has a strong influence. For this reason it is assumed in the lattice simulations that the bond elements have the same elastic modulus as the matrix, but a lower tensile strength. These assumptions were made on basis of results obtained from 2D lattice analyses using relatively coarse meshes (length of the lattice element 1 mm) and it is questionable whether the same assumptions hold when the ITZ is represented more realistically, by using finer lattice meshes. Here, the three cases already analysed in paragraph 5 are considered, where the bond strength is set equal to the matrix strength. The results are shown in Fig. 13. It is evident from the diagram which shows the number of beams failed during loading, that it is not possible anymore to distinguish the phase of de-bonding from the phase of crack propagation in the matrix. Rather, cracks form next to the grains or inside the matrix randomly. The effect on the load–displacement diagram is that the pre-peak region is nearly linear, does not contain any kinks. Consequently no variation of stiffness is observed. It should be mentioned that such a kink was earlier detected in 2D lattice analyses for a very small value of the bond strength [10]. Apparently, 3D effects seem to emphasize the non-linearity induced by de-bonding. Next, the particle density hardly affects the peak-load. Since the number of fractured elements decreases, the amount of dissipated energy decreases as well (in comparison to the case with low strength interface). Once again the load–displacements graphs can be better understood by considering the
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
939
Fig. 13. Results from lattice simulations with Pk;latt ¼ 29%, 39% and 53%, with ft;b ¼ ft;m : (a) force–displacement diagram; (b) cumulative number of beams removed versus displacement.
curves which show the number of beams removed at each loading stage. Although ITZ and matrix have exactly the same properties, they are still distinguished in these graphs, in order to identify their position with respect to the aggregates. The curves representing bond and matrix failure present a step approximately at the same loading stage, showing that the processes of de-bonding and crack formation and propagation through the matrix start and develop simultaneously. Furthermore, failure of the specimen occurs quicker than in the case considered in paragraph 5, i.e. less steps are necessary for reaching 25 lm displacement. Still, the number of matrix elements which fail decreases with increasing particle density, while the number of bond elements increases. Differently than in the case considered in paragraph 5, matrix failure prevails on bond failure. Only when Pk;latt ¼ 53% more de-bonding than matrix failure occurs. Even in this case, however, de-bonding is comparable to matrix failure. Also in the present case the fracture process is not confined to the matrix and bond phase, but also aggregate elements fail. The percentage of the total number of failed aggregate elements increases with the particle density, as in paragraph 5. However, the percentages are much higher now, namely 9%, 15% and 23% for Pk;latt ¼ 29%, 39%, and 53%, respectively (see Fig. 14(b)).
Fig. 14. Percentages of the total amount of failed beams for each phase: (a) the ITZ has a lower strength than the matrix; (b) the ITZ has the same strength as the matrix.
940
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
7. Conclusions In order to solve the problem of brittleness in force–displacement diagrams obtained from 2D lattice analyses a 3D lattice model was developed. Due to the enormous increase of number of degrees of freedom (and thus of computational time) the model was implemented in a finite element code for parallel computing. Numerical experiments were conducted on a cube subjected to uni-axial tension. The numerical results show that the model can reproduce the fracture processes observed in real physical experiments. Firstly, de-bonding is caused by stress concentrations around the particles where the material is weaker. Once the peak-load is reached cracks localize in the matrix. In the load–displacement diagram de-bonding corresponds to non-linear pre-peak behaviour, while crack formation and propagation in the matrix correspond to the softening branch. When the particle density increases, the fraction of bond elements increases and de-bonding prevails on matrix failure. As an effect on the load–displacement diagram, the peak-load decreases and the relative ductility increases. When the same strength of the matrix is assigned to the bond, de-bonding and matrix failure cannot be distinguished anymore. As a consequence, the pre-peak branch of the load–displacement diagram is nearly linear, the peak-load is hardly affected by the particle density. Furthermore, since the amount of dissipated energy is related to the number of fractured elements, it was found that, although increasing with particle density, the absolute value of the dissipated energy is smaller than in the case of low strength interface. Apparently, a 3D lattice version coupled with a more realistic representation of the particle distributions might be the solution to the problem of brittleness encountered in the 2D lattice analyses. However, due to limitations in time, the results presented in this paper are obtained each from a single lattice simulation. Since the model is statistical in nature, a more extensive series of numerical simulations would be necessary. At least, each analysis should be repeated three times, as is common practice in experiments. As a matter of fact, concrete is a material with a disordered internal structure. Because of this and possible boundary effects due to the relatively small size of the specimen considered, some scatter in the results must be obtained.
Acknowledgements The financial support by STW/PPM and the Ministry of Works (RWS-Bouwspeurwerk) is gratefully aknowledged, as is the expert help of Mr. Frank Everdij in performing the computations.
References [1] Hrennikoff A. Solution of problems of elasticity by the framework method. J Appl Mech (Trans ASME) 1941;8:A169–75. [2] Herrmann HJ, Hansen A, Roux S. Fracture of disordered, elastic lattices in two dimensions. Phys Rev B 1989;39(1):637–48. [3] Schorn H, Rode U. 3-D modelling of process zone in concrete by numerical simulation. In: Shah SP, Swartz SE, editors. Fracture of concrete and rock. New York: Springer Verlag; 1987. p. 220–8. [4] Wang J. Development and application of a micro-mechanics based numerical approach for the study of crack propagation in concrete. PhD thesis, Swiss Federal Institute of Technology, 1994. [5] Roelfstra PE, Sadouki H, Wittmann FH. Le beton numerique. Mater Struct (RILEM) 1985;18(107):327–35. [6] Vonk RA. Softening of concrete loaded in compression. PhD thesis, Eindhoven University of Technology, 1992. [7] Bolander Jr JE, Kobashi Y. Size effect mechanisms in numerical concrete fracture. In: Wittmann FH, editor. Fracture mechanics of concrete structures. AEDIFICATIO Publishers; 1995. p. 535–42. [8] Schlangen E, van Mier JGM. Experimental and numerical analysis of micro-mechanisms of fracture of cement-based composites. Cem Concr Compos 1992;14:105–18. [9] Lilliu G, van Mier JGM. Analysis of crack growth in the brazilian test. In: Eligehausen R, editor. Construction materials: theory and application (HW Reinhardt 60th birthday commemorative volume). Stuttgart: Ibidem Verlag; 1999. p. 123–37.
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
941
[10] Schlangen E. Experimental and numerical analysis of fracture processes in concrete. PhD thesis, Delft University of Technology, 1993. [11] Vervuurt A. Interface fracture in concrete. PhD thesis, Delft University of Technology, 1997. [12] van Vliet MRA. Size effect in tensile fracture of concrete and rock. PhD thesis, Delft University of Technology, 2000. [13] Lingen FJ. Design of an object oriented finite element package for parallel computers. PhD thesis, Delft University of Technology, 2000. [14] Moukarzel C, Herrmann HJ. A vectorizable random lattice. J Statist Phys 1992;68:911–23. [15] Lilliu G, van Mier JGM. Numerical characterization of the elastic properties of heterogeneous materials with a 3D lattice model. In: Carlomagno GM, Brebbia CA, editors. Computational methods and experimental measurements. Southampton (UK): WIT Press; 1999. p. 515–24. [16] Tzschichholz F. Peeling instability in cosserat like continua. Tech. Rep. Preprint HLRZ 68/91, HLRZ KFA J€ ulich, 1991. [17] van Mier JGM, Shi C. Stability issues in uniaxial tensile tests on brittle disordered materials. Int J Solids Struct 2002;39(13/ 14):3359–72. [18] Duan K, van Mier JGM. Crack growth in four concretes under monotonic or cyclic splitting load. In: Azizinamini A, Darwin D, French C, editors. Engineering foundation international conference on high strength concrete. Reston (VA): ASCE; 1999. p. 444– 56. [19] van Mier JGM, Chiaia B, Vervuurt A. Numerical simulation of chaotic and self-organizing damage in brittle disordered materials. Comput Meth Appl Mech Engng 1997;142(12):189–201. [20] Schlangen E, Garboczi EJ. New method for simulating fracture using an elastically uniform random geometry lattice. Int J Engng Sci 1996;34(10):1131–44. [21] Stroeven M. Discrete numerical modelling of composite materials. PhD thesis, Delft University of Technology, 1999. [22] Schlangen E, van Mier JGM. Lattice model for simulating fracture of concrete. In: Wittmann FH, editor. Numerical models in fracture mechanics of concrete. Rotterdam: Balkema; 1993. p. 195–205. [23] van Mier JGM. Crack face bridging in normal, high strength and lytag concrete. In: van Mier JGM, Rots JG, Bakker A, editors. Fracture processes in concrete rock and ceramics. London: Chapman & Hall/E&FN Spon; 1991. p. 27–40. [24] van Mier JGM. Mode I fracture of concrete: discontinuous crack growth and crack interface grain bridging. Cem Concr Res 1991;21(1):1–15.
3D lattice type fracture model for concrete G. Lilliu
a,*
, J.G.M. van Mier
a,b
a b
Microlab, Delft University of Technology, Delft, The Netherlands Institute for Building Materials, ETH, CH 8093 Z€urich, Switzerland
Received 21 February 2002; received in revised form 4 July 2002; accepted 17 July 2002
Abstract A 3D beam lattice model is used for simulating fracture processes in concrete, which is schematized as a three-phase material (aggregate, interfacial transition zone and matrix). Numerical experiments are conducted varying the particle density. The obtained results show that, with increasing particle density, the peak-load decreases and the response is more ductile. This appears to be related to the amount of de-bonding. As a matter of fact, when the strength of the interface transition zone is set equal to the strength of the matrix, neither the peak-load nor the ductility of the lattice response are influenced by the particle density. 2002 Elsevier Science Ltd. All rights reserved. Keywords: 3D lattice; Particle density; Peak-load; Ductility; Interfacial transition zone
1. Introduction Fracture processes due to mechanical loading are a major cause of damage in concrete structures. An adequate comprehension of these processes is important not only for ensuring safety of the construction, but also for optimizing the material behaviour, which would apport beneficial economical effects. Apparently, crack formation and propagation in concrete are affected by its highly disordered structure. Therefore, this structure should be taken into account modelling the material behaviour. At present fracture is modelled at different levels of observation, ranging from the macroscopic level, where material is considered homogeneous and non-linear constitutive relations are used, to meso, micro, or even to atomic or sub-atomic scale. At the meso-level of observation concrete is regarded as a two or three-phase material. A model which allows a straightforward implementation of the material heterogeneity at the meso-level is the lattice model. Such type of material modelling has already been used in the 1940s for simulating problems of elasticity [1]. More recently it has been introduced in the field of statistical physics for studying fracture of disordered materials [2]. Lattice models assume usually a linear-elastic material constitutive relation and, for this reason, they are conceptually very simple. Since the lattice models are a simple and straightforward tool for understanding the physics of fracture processes, they are considered particularly appealing to the authors.
*
Corresponding author. E-mail addresses: [email protected] (G. Lilliu), [email protected] (J.G.M. van Mier).
0013-7944/03/$ - see front matter 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 3 - 7 9 4 4 ( 0 2 ) 0 0 1 5 8 - 3
928
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
Such models, at least in their 2D version, present a characteristic brittleness of the response in terms of load–displacement diagrams. Apparently, this is a consequence of the assumed brittleness of the single element of the network which describes the material. However, it is also important to consider carefully to what extent the material is reproduced realistically in the model. In fact, in order to reduce the computational time, the adopted meshes are usually too coarse. The consequence of this is that only the largest fractions of the particle distributions are reproduced (with diameter da P 2 mm), while in real concrete grains with size of the order of micro-meters are present (e.g. cement grains). Furthermore, the adopted particle distributions, usually computer-generated, hardly are as dense as in real concrete. Last but not least, real fracture processes are 3D and not 2D. Presumably, both the presence of small particles and the 3D effect can increase the ductility of the lattice response. In fact, since more branching and bridging in the crack pattern is expected, the amount of energy dissipated in the post-peak region of the load–displacement diagram should increase. In the past many different types of meso-level models have been developed, for example [3–7]. In all these models the particle structure of concrete was included in a finite element mesh. In some of these models the interface between aggregate and matrix was modelled with zero-length spring elements. These elements had a normal softening component and a shear component along the surface of the aggregate [5,6]. The main objection against these models is that they are unsuitable to explain softening, since softening is already included at the meso-level. An alternative approach was that proposed by [4]. Based on linear-elastic fracture mechanics, crack propagation along the interface and the matrix could be described. The disadvantage of this model is that initial cracks must be defined at the beginning of the analysis, since no crack nucleation criterion is included. Lattice type models have been applied by others, too [3,7]. Schorn [3] used truss elements, which require some numerical measures to avoid instability during fracture propagation. Again, softening was introduced at element level in the model of Bolander [7], who applied a beam lattice with a step-wise softening law for the elements. In this paper results are shown of simulations made with a 3D version of the so-called ÔDelft latticeÕ. This is an extension of a previous 2D version [8], which has been extensively used in the past decade to simulate different experiments [9–12]. In this model beam elements are used for avoiding numerical instabilities. At element level, an elastic purely brittle fracture criterion is adopted, which makes the model suitable for studying the effect of material structure on softening. An additional advantage is that the particle structure can be included in a simple and straightforward manner. Since 3D lattice simulations with relatively fine meshes would be practically impossible with the standard software, the model has been implemented in a code for parallel computing [13]. The case considered is a prism subjected to uni-axial tension. The effect of varying particle density is studied for two cases: in the first case, the interfacial transition zone (ITZ) has a lower strength than the matrix; in the second case, the same mechanical properties as the matrix are assigned to the interface. Thus, in the latter case, the material is considered as two-phase material rather than containing three distinct phases.
2. The model The material is schematized with a network of Bernoulli beams. Lattices can be regular or random. In the first case, all elements of the mesh have the same length. In the second case, after the nodes have been positioned randomly inside the cells of a grid with regular cell size s, they are connected by means of the Voronoi construction [14] (see Fig. 1). For each node, a number of neighboring nodes is checked for possible connectivities (36 in the 2D case, 178 in the 3D case). In random lattices, where all elements of the mesh have different lengths and, thus, different stiffnesses, disorder is already introduced at the ‘‘geometrical’’ level. The amount of such disorder can be varied by placing the nodes inside sub-cells of the original cells. Randomness is defined as the ratio between the size of the sub-cell and the main cell, i.e. A=s. The cross-section of the beams is rectangular and circular in the 2D and 3D case, respectively. The height (in the
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
929
Fig. 1. (a) Voronoi construction for a 3 3 2D mesh with A=s ¼ 0:5; (b) generation of the connectivities for one node of the mesh shown in (c); 3 3 3 3D mesh with A=s ¼ 0:5.
2D case) or diameter (in the 3D case) of the cross-section, h, must be chosen properly in order to obtain the same overall PoissonÕs ratio as the real material [11,15]. Heterogeneities can be introduced in a straightforward manner. The mesh is overlaid on top of a computer-generated particle distribution and different mechanical properties are assigned to the beams falling in each phase. Commonly used ratios for the YoungÕs moduli and tensile strengths are Ea =Em ¼ 70 000=25 000, Em =Eb ¼ 25 000=25 000 and ft;a =ft;m ¼ 10=5, ft;m =ft;b ¼ 5=1:25 (where a, m and b refer to aggregate, matrix and bond, respectively). The resulting density of the aggregates (ratio between number of aggregate beams and total number of beams) is usually much lower than the particle density in real concrete. One of the reasons is that the particles are placed randomly, starting from the fraction with largest size, in a way that a minimum mutual distance is respected. As a consequence, placement of the smallest particle fractions becomes more and more difficult when the density increases. Another reason is that part of the particle volume falls inside the interface (bond zone), which is much thicker in the lattice than in reality. Finally, since the adopted meshes are relatively coarse, the smallest particles are automatically discarded (see Fig. 2(b)). In fact, the length of the beams and the minimum grain size, lbeam and da;min , should follow the minimum requirement lbeam 6 da;min =3. However, this would imply an enormous increase of the number of degrees of freedom and of the computational time. In lattice simulations, fracture processes are reproduced by sequential removal of elements from the mesh and it is assumed that such elements have linear-elastic behaviour until failure. At each step, the element to be removed is that with maximum value of the ratio between the effective stress r , and the tensile strength
Fig. 2. (a) Introduction of heterogeneities by overlaying the lattice mesh on top of a computer-generated particle distribution; (b) distinction between aggregate, matrix and bond phase according to the location of a lattice element in the particle structure.
930
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
Fig. 3. Smoothening of a force–displacement diagram obtained from a lattice simulation.
ft . Despite the assumed brittleness of local failure, global structural softening can be obtained. The effective stress is computed at the ends of the lattice beam, according to: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Mx2 þ My2 F ; r ¼ þ a A W where A and W are the area and the section modulus of the cross-section and F, Mx and My are the internal axial force and bending moments. Thus, it is assumed that only mode I fracture occurs locally. The coefficient a can be used to adjust the contribution of the bending moment in a beam. For 2D lattice analyses it was shown that small values of a give a long tail in the stress-deformation response [10]. However, if the value of a is taken equal to zero, peeling effects might occur in the fracture process in case of a rectangular lattice [16]. A value commonly used in 2D lattice simulations is a ¼ 0:005. It can be shown that in the 3D analyses reported in this paper a variation of a in the range 0.00–0.05 does not affect the crack pattern. Therefore, in this paper the value a ¼ 0 will be adopted, which means essentially that a normal stress criterion is used. For sake of simplicity, torsion is not included in the fracture law. Finally, a remark concerning the load–displacement diagrams resulting from lattice analyses should be made. They present a characteristic zig-zag pattern, which shows that local instabilities occur during crack formation and propagation. Nevertheless, the diagrams obtained from lattice simulations, with their drops
Fig. 4. Case study: (a) computer-generated particle distribution; (b) positioning of the specimen inside the volume Vtot .
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
931
of load, are quite similar to those recorded in the laboratory during experiments [17]. However, when a specific control method is adopted for conducting the experiment, stable crack propagation can be obtained. Still, local instabilities may occur, but they give rise to small indents in the curve only. The zig-zag in the lattice load–displacement diagrams is usually overcome by applying a smoothening technique as shown in Fig. 3. For comparison of various results, smoothening is quite helpful. The inconvenience of this procedure is that information about possible snap-backs is lost and that the amount of dissipated energy is overestimated. In the model the irreversible deformations due to geometrical mismatch of the crack faces during cyclic loading [18] are not taken into account. A way for including permanent deformations is to use a phenomenological approach like in [19].
3. Case study As case study a cube subjected to uni-axial tension is analysed (see Fig. 4(b)). The edge of the cube is 24 mm long. For generating the lattice mesh, the size s ¼ 0:75 mm of the gridpcells ffiffiffi was chosen, with randomness A=s ¼ 0. Therefore, the diagonal beams are the longest, with lbeam ¼ 3s ¼ 1:3 mm. This value of the randomness was selected in order to isolate the effect of the disorder induced by the particles from the geometrical disorder. The nodes of the upper and bottom face of the cube are supported in the X, Y, and Z direction, and an uniform displacement in the Y direction is applied to the nodes of the upper face. The choice of the geometry of the lattice was dictated only by the necessity to limit the computational time. In fact, the number of nodes in the mesh is 35 937, resulting in 431 244 degrees of freedom, which could be handled using 32 processors of a CRAY T3E. Obviously, such a small cube cannot be representative of the overall material behaviour, since boundary effects can be expected. Furthermore, only particles with da P 4 mm can be included, which are relatively big in comparison with the size of the specimen. In this case the material is not Ôstatistically homogeneousÕ and its internal structure (e.g. spatial distribution of particles) might play a major role, causing a large scatter of the outcome from the simulations. Nevertheless, such a numerical experiment is useful for a qualitative understanding of 3D fracture processes. An important parameter to be assigned in lattice analyses is the dimension of the cross-section, which influences the overall PoissonÕs ratio. In random lattices the value of m varies with the randomness of the mesh and when the randomness is A=s ¼ 0 never reaches the value 0.2 [11,15], which is the value usually adopted for concrete. However, it was found that m becomes independent from the randomness when h=s 0:8. For this reason, it was decided to adopt the same value h ¼ 0:577 mm like in simulations of 2D regular lattices with beam length lbeam ¼ 1 mm [20]. It should be mentioned that with the adopted geometry of the beams beam theory is not applicable any longer. This fact is ignored in the analyses to come. As a matter of fact the model is a very rough approximation of reality but this appears throughout the entire approach.
4. Fracture process with sparse particle distribution 4.1. Aggregate structure A particle distribution is considered with maximum grain size da;max ¼ 12 mm and volumetric density Pk ¼ 75%, which is realistic for concrete. A Fuller grading curve is adopted as basis for the computergenerated aggregate structure. The number of particles with diameter di1 6 da 6 di is obtained following qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi di =da;max di1 =da;max Vtot ; Pk
932
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
Fig. 5. (a) Grading curves of the computer-generated particle distribution shown in Fig. 4(a); (b) phase fractions after overlaying the particle structure on top of the 3D lattice mesh with randomly varying position of the specimen inside Vtot .
where Vtot is the volume to be filled (in the example Vtot ¼ 72 72 72 mm3 ). In order to discretize the Fuller curve, the range 0 mm 6 da 6 12 mm of diameters is divided into intervals of 0.5 mm, as shown in Fig. 5(a). As an approximation, all grains falling in one interval of the distribution are considered to have the same diameter, namely, the integer value of the interval. As a consequence, the grains falling in the interval 0 mm 6 da 6 0:5 mm are neglected and the resulting volumetric density decreases to Pk;eff ¼ 60%. Next, the grains are positioned randomly inside the volume Vtot . Starting from the biggest diameter, the particles are placed in a way that a mutual minimum distance, 1:1ðda;i þ da;j Þ=2, between two adjacent particles is guaranteed. Since positioning of particles becomes difficult when the smallest fractions are considered they are usually neglected. Finally, geometrical considerations (see paragraph 2) impose a lower bound to the minimum diameter of the particles (4 mm 6 da ), causing a further decrease of density, until the value Pk;eff ¼ 34%. Once the particle distribution has been generated, the specimen can be placed in the volume Vtot (see Figs. 2(a) and 4 for a 2D and a 3D case, respectively) and heterogeneities can be implemented. After the lattice has been overlaid on top of the particles a new density Pk;latt can be defined as the ratio between the number of aggregate beams and the total number of beams. Fig. 5(b) shows how the fractions of the three phases vary with randomly varying the position of the specimen. As already men tioned in paragraph 2, Pk;latt 6 Pk;eff . In paragraph 4.2 the results obtained from a case where Pk;latt ¼ 29% will be reported. 4.2. Results For the case Pk;latt ¼ 29%, the smoothened force–displacement diagram, until a displacement of 25 lm is shown in Fig. 6(a). A detail of the same diagram, until 5 lm displacement, is shown in the inset. The diagram shows two main branches, namely, the pre-peak and the post-peak branch. In the pre-peak regime, after the first beam element has failed (at the stage marked with a, which corresponds to 20% of the peakload), the diagram shows a kink, which is followed by a nearly linear part (b–c). A strong non-linearity (c– d) can be observed before the peak-load. In the post-peak regime the diagram shows subsequently a steep drop of load (until the stage marked with e) and a long tail. In order to enlighten the development of the fracture process the number of beams belonging to matrix or bond phase, which have failed at each loading stage, is shown in Fig. 6(b). Also aggregate beams fail during the fracture process. However, only 2% of the total number of beams failed at 25 lm displacement are aggregate elements. This value is negligible when compared with the percentages of failed matrix and bond elements, which are 27% and 71%, respectively
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
933
Fig. 6. (a) Force–displacement diagram obtained from a lattice simulation with Pk;latt ¼ 29%; (b) cumulative number of failed beams versus displacement.
(see Fig. 14(a)). Firstly, it can be observed that more interface elements rather than matrix elements fail. Although the first matrix beam fails at point b (60% of the peak-load), practically only de-bonding occurs until point c of the diagram is reached (85% of the peak-load). As a matter of fact, the reason for debonding to prevail is twofold: • stress concentrations occur next to the particles due to the difference in stiffness between the aggregate and the bond/matrix phases, and, • the interface has the lowest strength in the three-phase material. After point c, cracks start to appear in the matrix and propagate rapidly after the peak-load. This corresponds to a steep drop in load, directly following the maximum load. The same phenomenon can be observed in the laboratory and corresponds to macroscopic crack propagation. Finally, prevailing failure of matrix elements is detected along the tail of the softening branch of the diagram. For visualizing the crack process, the crack patterns at stages d, e and at the value 25 lm are shown in Fig. 7. There, an exploded view of the specimen, after it has been cut into 0.75 mm thick slices, is given. The failed beams, falling inside each slice, are projected on a plane, and represented in white or black, depending if they are interface or matrix beams, respectively. From Fig. 7 three stages can be clearly distinguished: • first stage: cracking is due mostly to distributed de-bonding, • second stage: major cracks have formed, • third stage: no relevant changes in the fracture process occur anymore, but crack localization becomes more pronounced.
5. Effect of the particle density on the fracture process 5.1. Aggregate structure As already mentioned in the previous paragraph, high packing densities cannot be obtained with random positioning of the grains (unless much time consuming methods are used). Although yet not as high as in reality, the density of the computer-generated particle distributions can be increased reducing gradually the volume containing the initial sparse packing, by means of dynamic impact models. Using the program
934
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
Fig. 7. Exploded view of the crack patterns in the specimen with Pk;latt ¼ 29% at different loading stages: (a) at peak-load; (b) at 4 lm displacement (point e in the force–displacement diagram); (c) at 25 lm displacement (end of the tail in Fig. 6(a)). The failed bond and matrix elements are represented in white and black, respectively.
SPACE [21], which follows such an approach, it was possible to generate particle distributions with Pk;eff ¼ 48% and 62%, respectively. The difference with the particle distribution mentioned in paragraph 4 is that here the diameter of the particles varies continuously in the range 4 mm 6 da 6 12 mm. Fig. 8(a) shows the number of particles belonging to each fraction of the sieve curve, for the three different densities. In Fig. 8(b) the variation of the phase fractions with the density are shown. For computing the number of beams falling in each phase, one position of the specimen inside the volume containing the particles has been adopted (the same position like that considered in the previous paragraph). Once the lattice is overlaid on
Fig. 8. (a) Number of particles for three different values of the volumetric aggregate density, Pk;eff ; (b) phase fractions of the 3D lattice . as function of Pk;eff
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
935
top of the particle distribution, the densities Pk;latt ¼ 39% and 53%, corresponding respectively to Pk;eff ¼ 48% and 62%, are obtained. The matrix phase decreases with increasing particle density, while the bond phase increases, although slower than the aggregate fraction. The two densities, Pk;latt ¼ 39% and 53%, correspond to the two extreme situations of comparable amount of aggregate and matrix phase, and comparable amount of matrix and bond phase, respectively. 5.2. Results The outcome from the lattice simulations for the values of density Pk;latt ¼ 29%, 39% and 53%, is shown in the diagrams of Fig. 9. The dimensionless load–displacement diagrams, where load and displacement are divided by the values at the peak, are shown in the inset of Fig. 9(a). The number of failed beams versus displacement is shown in Fig. 9(b). Finally, the crack patterns for Pk;latt ¼ 39% and 53%, respectively, are shown in Figs. 10 and 11. From such figures it appears that both the peak-load and the post-peak behaviour are influenced from the particle density. Nevertheless, the fracture process still develops as described in the previous paragraph. Initially, only de-bonding occurs. Next, cracking localizes in the matrix (see Fig. 9(b)) and longer cracks develop in the specimen (see Figs. 7, 10 and 11). Crack localization corresponds to the initial part of the post-peak branch in the load–displacement diagram. A relevant difference among the three cases should be mentioned. For Pk;latt ¼ 29% and 39%, crack localization corresponds to a steep branch after the peak, and matrix failure prevails on interface failure in the tail of the softening curve. For Pk;latt ¼ 53%, after the peak the load decreases gradually while the displacement increases, which produces a more ductile response. Furthermore, interface failure still prevails in the tail of the softening curve. In the next two paragraphs, the effect of the particle density on the peak-load and the ductility of the lattice response, respectively, will be discussed in more detail. 5.2.1. Effect on the peak-load It can be observed from Fig. 9(a) that the peak-load decreases with particle density. This can be explained as follows. The carrying capacity of the specimen depends on the strength of the matrix–bond system. The interface fraction is the weakest component of such a system and contributes to lower its strength. Thus, since the fraction of bond elements in the lattice increases with the density, it can be expected that lower peak-loads are achieved with higher particle densities. Also another factor can contribute
Fig. 9. (a) Force–displacement diagram obtained from lattice simulations with Pk;latt ¼ 29%, 39% and 53%; (b) cumulative number of removed beams versus displacement.
936
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
Fig. 10. Exploded view of the crack patterns in the specimen with Pk;latt ¼ 39% at different loading stages: (a) at peak-load; (b) at 4 lm displacement; (c) at 25 lm displacement.
Fig. 11. Exploded view of the crack patterns in the specimen with Pk;latt ¼ 53% at different loading stages: (a) at peak-load; (b) at 4 lm displacement; (c) at 25 lm displacement.
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
937
to this. At the peak-load the fracture process starts to localize in the matrix. At this stage the structure contains a certain amount of damage, which is related to the amount of de-bonding in the pre-peak regime. As a matter of fact, with increasing number of bond elements more of such elements are removed before cracks start to propagate in the matrix (see Fig. 9(b)). The consequence of this is that the structure is more damaged and its carrying capacity decreases. The dependence of the peak-load on the volume and strength of bond phase will be discussed again in paragraph 6. There it will be shown that, if no bond phase is considered, the peak-load does not depend on the particle density. 5.2.2. Effect on the ductility Here, ductility is defined as the ratio between the energy consumed during the fracture process after and before the peak. The results show that an increase of the particle density produces an up-lift of the softening branch in the dimensionless load–displacement diagram (see inset in Fig. 9) and, thus, an increase of ductility. This is more remarkable for the highest value of the density, Pk;latt ¼ 53%, which approaches the real particle density in concrete (75%). In this case de-bonding plays an important role also in the postpeak region. In this region, not only the number of failing bond elements increases continuously, but even with a higher rate than the matrix elements (see Fig. 9(b)). As pointed out in the introduction, a higher ductility is expected when more branching and bridging occur [22]. Both these phenomena were observed in impregnation experiments during the tail of the softening curve [23,24]. Indeed the current analyses show that more branching, which implies more de-bonding, contribute to increase the ductility of the lattice response. Due to the difficulty in visualizing the crack patterns resulting from the simulations, it is harder to verify that also bridging can be reproduced in lattice simulations. However, matrix failure in the postpeak regime might be related to flexural cracking due to bridging, and not only of macroscopic crack propagation. Thus, realistic particle distributions seem to be essential to obtain a more ductile response. However, some comments should be made. In Fig. 14, the percentage of failed elements falling inside each of the three phases is represented as function of the density (such percentage is computed on the total number of elements removed at a displacement of 25 lm). Fig. 14(a) shows that the percentage of aggregate elements failing during the experiment increases with the density. For Pk;latt ¼ 53%, such percentage reaches the value 4.6%, while the percentages of failed matrix and bond elements are 13.7% and 81.7%, respectively. Although the percentage of failed aggregate beams is negligible when compared with the percentage of failed interface elements, it cannot be ignored when compared with the percentage of failed matrix elements. It is questionable whether the fracture of aggregate beams can affect somehow the ductility of the lattice response. It might be expected that in normal concrete, like is the material considered here, fracture is limited only to interface and matrix. Nevertheless, since also in experiments aggregate failure sometimes occurs, the lattice results still appear quite realistic. The fact that also aggregate elements are fractured might be related to stress concentrations induced by the geometry of the grains. As a matter of fact, since the lattice meshes are always relatively coarse, the contour of the grains is never smooth but shows irregularities (see Fig. 12(b)). Most probably, a more appropriate choice of mechanical properties would have avoided aggregate fracture. For example, it might have been sufficient to assume a higher aggregate strength, i.e. ft;a ¼ 4ft;m rather than ft;a ¼ 2ft;m .
6. Effect of the bond strength on the fracture process From the results obtained with varying particle density it appears that the ITZ, which is the weak layer linking the aggregates to the matrix, plays a major role in the fracture process. The importance of the ITZ is evident also from Fig. 12(b), where a detail of a 2D regular lattice with beam length 0.25 mm is shown. In Fig. 12(b) the particle distribution has a density Pk;eff ¼ 83%, with minimum diameter 1 mm and maximum diameter 12 mm. For such a high density, even though a very fine lattice mesh is used, continuous paths of
938
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
Fig. 12. (a) Schematization of the ITZ in the material with a series model; (b) detail from a lattice mesh where percolation of the bond phase occurs.
bond elements can be detected. In such case the cracks in the bond phase will coalesce, causing formation of macrocracks (percolation). There is quite some debate on the thickness of the real ITZ, which varies in the range between 30 and 50 lm. In the current model the bond zone is a relatively wide zone (at least 10 times as large as the real ITZ), which should be better regarded as a region of average properties among the three phases, namely aggregate, ITZ and matrix. The bond element can be schematized in a simple manner adopting a system of longitudinal springs connected in series (see Fig. 12(a)), with lengths la , lb , and lm , respectively. For such a system, the equivalent elastic modulus Eb;eq , follows l la lb lm ¼ þ þ Eb;eq Ea Eb Em and the equivalent strength is the same as the ITZ strength which is the weakest element of the system. In [11] it was shown that the elastic modulus of the ITZ hardly affects the fracture process, while the strength has a strong influence. For this reason it is assumed in the lattice simulations that the bond elements have the same elastic modulus as the matrix, but a lower tensile strength. These assumptions were made on basis of results obtained from 2D lattice analyses using relatively coarse meshes (length of the lattice element 1 mm) and it is questionable whether the same assumptions hold when the ITZ is represented more realistically, by using finer lattice meshes. Here, the three cases already analysed in paragraph 5 are considered, where the bond strength is set equal to the matrix strength. The results are shown in Fig. 13. It is evident from the diagram which shows the number of beams failed during loading, that it is not possible anymore to distinguish the phase of de-bonding from the phase of crack propagation in the matrix. Rather, cracks form next to the grains or inside the matrix randomly. The effect on the load–displacement diagram is that the pre-peak region is nearly linear, does not contain any kinks. Consequently no variation of stiffness is observed. It should be mentioned that such a kink was earlier detected in 2D lattice analyses for a very small value of the bond strength [10]. Apparently, 3D effects seem to emphasize the non-linearity induced by de-bonding. Next, the particle density hardly affects the peak-load. Since the number of fractured elements decreases, the amount of dissipated energy decreases as well (in comparison to the case with low strength interface). Once again the load–displacements graphs can be better understood by considering the
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
939
Fig. 13. Results from lattice simulations with Pk;latt ¼ 29%, 39% and 53%, with ft;b ¼ ft;m : (a) force–displacement diagram; (b) cumulative number of beams removed versus displacement.
curves which show the number of beams removed at each loading stage. Although ITZ and matrix have exactly the same properties, they are still distinguished in these graphs, in order to identify their position with respect to the aggregates. The curves representing bond and matrix failure present a step approximately at the same loading stage, showing that the processes of de-bonding and crack formation and propagation through the matrix start and develop simultaneously. Furthermore, failure of the specimen occurs quicker than in the case considered in paragraph 5, i.e. less steps are necessary for reaching 25 lm displacement. Still, the number of matrix elements which fail decreases with increasing particle density, while the number of bond elements increases. Differently than in the case considered in paragraph 5, matrix failure prevails on bond failure. Only when Pk;latt ¼ 53% more de-bonding than matrix failure occurs. Even in this case, however, de-bonding is comparable to matrix failure. Also in the present case the fracture process is not confined to the matrix and bond phase, but also aggregate elements fail. The percentage of the total number of failed aggregate elements increases with the particle density, as in paragraph 5. However, the percentages are much higher now, namely 9%, 15% and 23% for Pk;latt ¼ 29%, 39%, and 53%, respectively (see Fig. 14(b)).
Fig. 14. Percentages of the total amount of failed beams for each phase: (a) the ITZ has a lower strength than the matrix; (b) the ITZ has the same strength as the matrix.
940
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
7. Conclusions In order to solve the problem of brittleness in force–displacement diagrams obtained from 2D lattice analyses a 3D lattice model was developed. Due to the enormous increase of number of degrees of freedom (and thus of computational time) the model was implemented in a finite element code for parallel computing. Numerical experiments were conducted on a cube subjected to uni-axial tension. The numerical results show that the model can reproduce the fracture processes observed in real physical experiments. Firstly, de-bonding is caused by stress concentrations around the particles where the material is weaker. Once the peak-load is reached cracks localize in the matrix. In the load–displacement diagram de-bonding corresponds to non-linear pre-peak behaviour, while crack formation and propagation in the matrix correspond to the softening branch. When the particle density increases, the fraction of bond elements increases and de-bonding prevails on matrix failure. As an effect on the load–displacement diagram, the peak-load decreases and the relative ductility increases. When the same strength of the matrix is assigned to the bond, de-bonding and matrix failure cannot be distinguished anymore. As a consequence, the pre-peak branch of the load–displacement diagram is nearly linear, the peak-load is hardly affected by the particle density. Furthermore, since the amount of dissipated energy is related to the number of fractured elements, it was found that, although increasing with particle density, the absolute value of the dissipated energy is smaller than in the case of low strength interface. Apparently, a 3D lattice version coupled with a more realistic representation of the particle distributions might be the solution to the problem of brittleness encountered in the 2D lattice analyses. However, due to limitations in time, the results presented in this paper are obtained each from a single lattice simulation. Since the model is statistical in nature, a more extensive series of numerical simulations would be necessary. At least, each analysis should be repeated three times, as is common practice in experiments. As a matter of fact, concrete is a material with a disordered internal structure. Because of this and possible boundary effects due to the relatively small size of the specimen considered, some scatter in the results must be obtained.
Acknowledgements The financial support by STW/PPM and the Ministry of Works (RWS-Bouwspeurwerk) is gratefully aknowledged, as is the expert help of Mr. Frank Everdij in performing the computations.
References [1] Hrennikoff A. Solution of problems of elasticity by the framework method. J Appl Mech (Trans ASME) 1941;8:A169–75. [2] Herrmann HJ, Hansen A, Roux S. Fracture of disordered, elastic lattices in two dimensions. Phys Rev B 1989;39(1):637–48. [3] Schorn H, Rode U. 3-D modelling of process zone in concrete by numerical simulation. In: Shah SP, Swartz SE, editors. Fracture of concrete and rock. New York: Springer Verlag; 1987. p. 220–8. [4] Wang J. Development and application of a micro-mechanics based numerical approach for the study of crack propagation in concrete. PhD thesis, Swiss Federal Institute of Technology, 1994. [5] Roelfstra PE, Sadouki H, Wittmann FH. Le beton numerique. Mater Struct (RILEM) 1985;18(107):327–35. [6] Vonk RA. Softening of concrete loaded in compression. PhD thesis, Eindhoven University of Technology, 1992. [7] Bolander Jr JE, Kobashi Y. Size effect mechanisms in numerical concrete fracture. In: Wittmann FH, editor. Fracture mechanics of concrete structures. AEDIFICATIO Publishers; 1995. p. 535–42. [8] Schlangen E, van Mier JGM. Experimental and numerical analysis of micro-mechanisms of fracture of cement-based composites. Cem Concr Compos 1992;14:105–18. [9] Lilliu G, van Mier JGM. Analysis of crack growth in the brazilian test. In: Eligehausen R, editor. Construction materials: theory and application (HW Reinhardt 60th birthday commemorative volume). Stuttgart: Ibidem Verlag; 1999. p. 123–37.
G. Lilliu, J.G.M. van Mier / Engineering Fracture Mechanics 70 (2003) 927–941
941
[10] Schlangen E. Experimental and numerical analysis of fracture processes in concrete. PhD thesis, Delft University of Technology, 1993. [11] Vervuurt A. Interface fracture in concrete. PhD thesis, Delft University of Technology, 1997. [12] van Vliet MRA. Size effect in tensile fracture of concrete and rock. PhD thesis, Delft University of Technology, 2000. [13] Lingen FJ. Design of an object oriented finite element package for parallel computers. PhD thesis, Delft University of Technology, 2000. [14] Moukarzel C, Herrmann HJ. A vectorizable random lattice. J Statist Phys 1992;68:911–23. [15] Lilliu G, van Mier JGM. Numerical characterization of the elastic properties of heterogeneous materials with a 3D lattice model. In: Carlomagno GM, Brebbia CA, editors. Computational methods and experimental measurements. Southampton (UK): WIT Press; 1999. p. 515–24. [16] Tzschichholz F. Peeling instability in cosserat like continua. Tech. Rep. Preprint HLRZ 68/91, HLRZ KFA J€ ulich, 1991. [17] van Mier JGM, Shi C. Stability issues in uniaxial tensile tests on brittle disordered materials. Int J Solids Struct 2002;39(13/ 14):3359–72. [18] Duan K, van Mier JGM. Crack growth in four concretes under monotonic or cyclic splitting load. In: Azizinamini A, Darwin D, French C, editors. Engineering foundation international conference on high strength concrete. Reston (VA): ASCE; 1999. p. 444– 56. [19] van Mier JGM, Chiaia B, Vervuurt A. Numerical simulation of chaotic and self-organizing damage in brittle disordered materials. Comput Meth Appl Mech Engng 1997;142(12):189–201. [20] Schlangen E, Garboczi EJ. New method for simulating fracture using an elastically uniform random geometry lattice. Int J Engng Sci 1996;34(10):1131–44. [21] Stroeven M. Discrete numerical modelling of composite materials. PhD thesis, Delft University of Technology, 1999. [22] Schlangen E, van Mier JGM. Lattice model for simulating fracture of concrete. In: Wittmann FH, editor. Numerical models in fracture mechanics of concrete. Rotterdam: Balkema; 1993. p. 195–205. [23] van Mier JGM. Crack face bridging in normal, high strength and lytag concrete. In: van Mier JGM, Rots JG, Bakker A, editors. Fracture processes in concrete rock and ceramics. London: Chapman & Hall/E&FN Spon; 1991. p. 27–40. [24] van Mier JGM. Mode I fracture of concrete: discontinuous crack growth and crack interface grain bridging. Cem Concr Res 1991;21(1):1–15.