Fracturing in concrete via lattice-particle model

II International Conference on Particle-based Methods - Fundamentals and Applications PARTICLES 2011 E. O˜nate and D.R.J. Owen (Eds) FRACTURING IN CONCRETE VIA LATTICE-PARTICLE MODEL JAN ELI´ Aˇ S∗AND ZDENˇ EKP.BA ˇ ZANT† ∗Masaryk-Fulbright Fellow, Northwestern University; Assistant Professor on leave from the Brno University of Technology, Czech Republic. e-mail: [email protected] †McCormick School Professor and W.P. Murphy Professor of Civil Engineering and Materials Science, Northwestern University, 2145 Sheridan Rd., CEE, Evanston, Illinois 60208; corresponding author, e-mail: [email protected] Key words: fracture, concrete, lattice-particle model, size effect, notch variability Abstract. Numerical simulation is used to explore the behavior of concrete beams of different sizes and different notch lengths, loaded in three-point bending. The entire range of notch depth is studied. One limit case is type 1 fracture, which occurs when the notch depth is zero and the crack initiates from a smooth surface (this is the case of the modulus of rupture test). Another limit is type 2 fracture, which occurs for deep enough notches. Both cases exhibit very different size effects. The fracture is simulated numerically with a robust mesolevel lattice-particle model. The results shed light on the transitional behavior in which the notch depth is non-zero but not deep enough for developing the the type 2 size effect dominated by energy release from the structure. In agreement with experimental observations and theoretical predictions, the numerical results show evidence of a decreasing macroscopic fracture energy as the ligament gets very short. 1 INTRODUCTION Modeling of the initiation and propagation of cracks in quasibrittle materials exhibiting strain softening has been studied for several decades. Although this is a difficult task complicated by the distributed damage dissipating energy within a fracture process zone (FPZ) of non-negligible size, realistic results have been achieved by some approaches; see e.g. [1]. In this contribution, the fracturing in concrete is modelled by the lattice-particle model developed in [2, 3, 4]. The main goal is to describe the transition between two basic types of failure. In type 1, a macroscopic crack initiates from a smooth surface and, in type 2, the crack 1 306 Jan Eli´aˇs and Zdenˇek P. Baˇzant initiates from a sufficiently deep notch or preexisting fatigued (stress-free) crack. Simple laws giving good approximations of test data have been derived for both types. However, the transition between these two types in the case of very shallow notches remains to be a challenge. In type 1, a large zone of distributed fracturing develops at the smooth surface until the damage localizes into a crack in the statistically weakest place described by the weakestlink model. In type 2, by contrast, the location of the crack is not random and a much smaller zone of distributed damage grows at a fixed place, the notch tip, until a state of critical energy release rate is reached. As the notch is getting shallower, the size of the damage zone increases and the crack location gradually develops random scatter. However, up to now there exists no experimental evidence for this transition, and so experiments to characterize it are in preparation at Northwestern University. The present purpose is to clarify this transition by numerical simulations, considering geometrically similar three-point bend concrete beams of constant thickness b, various depths Dand various relative notch depths α0. The present analysis is based on the cohesive crack model [5, 6, 1] (also called fictitious crack model). In this model, it is assumed that the cohesive stress transmitted across the crack is released gradually as a decreasing function of the crack opening, called the cohesive softening curve. Its main characteristic is the total fracture energy, GF– a material constant representing the area under this curve. For stationary propagation, the J-integral shows that GFalso represents the flux of energy into the FPZ. The fracture energy dissipation occurs within numerous meso-level microcracks in the FPZ. The present numerical model will directly simulate the behavior of these microcracks on the meso-level of a brittle inhomogeneous material such as concrete. For this purpose, the present analysis will be based on the discrete lattice-particle developed by G. Cusatis and coworkers [4], which is an extension of [2, 3]. The meso-level material fracture properties are in this model characterized by stress-displacement relations at the interfaces between grains or particles, representing the mineral aggregates in concrete. 2 BRIEF MODEL DESCRIPTION The material is represented as a discrete three-dimensional assembly of rigid cells. The cells are created by tessellation according to pseudo-random locations and radii of computer generated grains/particles. Every cell contains one grain (Fig. 1a,b). On the level of rigid cell connection, the cohesive crack model is used to represent the cracking in the matrix between the adjacent grains. The fracture energy is the same for all connections except that it depends on the direction of straining. The inter-particle fracturing is assumed to be of damage mechanics type. Thus the plastic frictional slip is not separately accounted for. But this simplification would matter only for unloading behavior which is not the objective of the present analysis. For a detailed description of the connection constitutive law or other model features, see [4]. However, the following minor deviations from the model in [4] are introduced: 2 307 Jan Eli´aˇs and Zdenˇek P. Baˇzant Figure 1: a) A rigid cell surrounding a particle and b) its section through the particle center; c) specimen geometry; d) coupling of the lattice-particle model with the standard elastic finite elements; e) example of pre-peak and post-peak damage pattern for the modulus of rupture test (zero notch depth). •The interparticle connection cohesive law in tension and shear is bilinear instead of exponential. It is therefore defined by eight constants: i) initial mesolevel fracture energy in tension and pure shear, Gtand Gs; ii) total mesolevel fracture energy in tension and pure shear, GTand GS; iii) the mesolevel cohesive tensile and shear strengths, σtand σs, and iv) the coordinates of the ”knee point”, i.e., the intersection of the two linear segments considered as 20% of tensile strength σtor shear strength σs, respectively. •The notch is represented simply by removing all the connections that cross the midspan provided that at least one of the centers of the connected particles is closer to the crack mouth than α0D. The advantage of this approach is that all grain positions can be completely random, and that the cutting of the notch by a saw is represented faithfully. The disadvantages are that the notch tip location is not exact and it is impossible to introduce notches whose depth is less than the minimal grain radius. •The confinement effect is neglected, but it was estimated that, in this type of experiment, the confinement does not play any important role. The mesolevel material properties in this model are deterministic. Randomness is introduced solely by pseudo-random locations and radii of grains. The effect of spatial 3 308 Jan Eli´aˇs and Zdenˇek P. Baˇzant variability of the material properties, which was found to be very important for capturing the statistical (Weibull) part of the type 1 size effect [7, 8]), is neglected. Since all the interparticle connections have identical deterministic fracture energy and tensile strength, the crack initiation from a smooth surface is preceded by distributed fracturing along the entire bottom surface. Nevertheless, the localized macroscopic crack always initiates very close to the midspan (Fig. 1e). 3 SIMULATION OF BEAMS OF VARIOUS SIZES AND NOTCH DEPTHS Beams geometrically similar in two dimensions, having depths D= 100, 200, 300, 400 and 500 mm and the same thickness of t= 0.04 m, were modelled. The span-depth ratio was S/D =2.4, and the maximal aggregate diameter was 9.5 mm. The minimal grain diameter was chosen as 3 mm. Based on the Fuller curve, particles of radii within chosen range were generated and pseudo-randomly placed into the specimen domain. The parameters of the connection constitutive law, which were mostly taken similar to those in [4], were: Ec= 30 GPa; Ea= 90 GPa; σt=2.7 MPa; Gt= 15 N/m; GT= 30 N/m; σs=3σt=8.1 MPa; Gs= 215 N/m; GS= 430 N/m; σc= 16σt= 43.2 MPa; Kc= 7.8 GPa; α=0.15; β= 1; µ=0.2; nc= 2. To ensure numerical stability in presence of softening, the simulations were controlled by prescribing the increase of the crack mouth opening displacement (CMOD) in every step. For unnotched beams, the location of macrocrack initiation was not known in advance, and so the controlling displacement was chosen to span several maximum aggregate sizes along the tensile face of the beam. To save computer time, the lattice-particle model covered only the region in which cracking was deemed to be possible. The region in which no damage was expected was assumed to follow linear elasticity and was modelled by standard 8-node isoparametric finite elements. The elastic constants for these elements were identified by fitting a displacement field with homogeneous strain to the discrete field of particle displacements generated at low stress level for a prism of particles subjected to low-level uniaxial compression. The macroscopic Young’s modulus and Poisson ratio were thus found to be E= 30.3 GPa and ¯ν=0.225. The finite element mesh was connected to the system of particles by introducing interface nodes treated as auxiliary zero-diameter particles (Fig. 1d). Same as the standard particles of the lattice model, these auxiliary particles had three translational and three rotational degrees of freedom. Each auxiliary particle lied at the boundary of one finite element. A similar interfacing was used in [3] but here, in contrast, the FEM nodes were considered to be the masters, and the auxiliary particle displacements were dictated by the master displacements according to the master element shape (or interpolation) functions. The rotations of the auxiliary particles were unconstrained. For large specimens and shallow notches, many particles are needed to fill the damage region. This led to extreme computational time and memory requirements. Therefore, such simulations were terminated as soon as the load dropped to 90% of the peak force. 4 309 Jan Eli´aˇs and Zdenˇek P. Baˇzant Figure 2: Damage patterns at the maximum load for beams of size D= 300 mm. Thin straight lines inside specimens show interface between lattice-particle model and finite elements. The benefit was that there was no need to simulate the crack in the upper part of the beam, and thus many fewer particles were needed. For unnotched beams or beams with deep notches, respectively, the simulations were run until the load was reduced to 40% or 5% of the peak load. Six realizations were computed for unnotched specimens of each size and for relative notch depth α0=a/D =0.5, but only three realizations for other notch depths, i.e., for α0∈{0.01, 0.02, 0.03, 0.05, 0.1, 0.2, 0.3, 0.4, 0.6, 0.7, 0.8, 0.9}. According to a procedure developed in [4], the depth of specimen was divided into horizontal strips of width h= 12 mm (Fig. 1c). The average stress σxnormal to the crack plane and the energy gddissipated per unit area were measured and stored for each strip and each step. This allowed reconstructing the macroscopic cohesive softening law in postprocessing. The σxvalue was obtained as the stress that must be applied on a vertical section through the center of each cell (Fig. 1b) to ensure equilibrium of both parts [9]. 4 OBSERVED NOTCHED-UNNOTCHED TRANSITION The simulations of beams with no notch (type 1) and deep notch (type 2) behaved as expected. A large fracturing zone developed before the peak load in the case of initiation from a smooth surface (type 1) (Fig. 2). For deep notches, the damage was localized above the notch tip only. For shallow notches, the damage above the tip prevailed but some mesolevel cracking along the bottom surface took place as well. A qualitatively different damage pattern at peak load was obtained comparing the shallowest notch and no notch simulations (Fig. 2 α0= 0 and 0.01). The smallest size D= 100 mm and smallest notch α0=0.01 (for which the actual notch depth is only 1 mm) was an exception. The procedure used for notch creation did not lead to the removal of any contact and so no notch was actually introduced. Similarly, a notch of 5 310 Jan Eli´aˇs and Zdenˇek P. Baˇzant 0 20 40 60 80 0 50 100 150 200 250 300 350 400 450 500 distance from bottom surface [mm] dissipated energy [N/m] a =0.5 0 a =0.7 0 a =0.3 0 a =0.2 0 a =0.1 0 a =0.05 0 a =0 0 a =0.01 0 0.02 0 20 40 60 80 0 50 100 150 200 250 300 350 400 distance from bottom surface [mm] dissipated energy [N/m] a =0.5 0 a =0.7 0 a =0.3 0 a =0.2 0 a =0.1 0 a =0.05 0 a =0 0 a =0.01 0 0.02 0 20 40 60 80 0 50 100 150 200 250 300 distance from bottom surface [mm] dissipated energy [N/m] a =0.5 0 a =0.7 0 a =0.3 0 a =0.2 0 a =0.1 0 0.05 a =0 0 a =0.01 0 0.02 0 50 100 150 200 a =0.5 0 a =0.7 0 a =0.3 0 a =0.2 0 a =0 0 a =0.01 0 0.1 0.05 0.02 0 20 40 60 80 0 50 100 dissipated energy [N/m] 0.3 a =0.05 0 0.1 0.2 a =0-0.02 0 0 20 40 60 80 dissipated energy [N/m] D = 500 mm D = 400 mm D = 300 mm D = 200 mm D = 100 mm Figure 3: Sums of energy per unit area dissipated in strips until stress σxin each strip gets reduced in post-peak to 75% of the tensile strength σt. the depth of 2 mm (obtained for D= 100 mm, α0=0.02 and D= 200 mm, α0=0.01) was represented poorly because only a few contacts, discontinuous through the specimen width, were severed. The simulation with mesolevel refinement of the microstructure does not suffice for capturing the effect of notches smaller than the dominant mesolevel inhomogeneities. A deeper refinement of microstructure simulation would be needed to capture the precise notch tip, but it would considerably complicate the programming. Based on the observed damage patterns at the peak, the type 1–type 2 transition begins immediately upon introducing any notch exceeding the size of the dominant inhomogeneity (>3 mm). Considering that the pure type 2 behavior is reached when almost no cracking occurs at the bottom surface, the transition becomes complete when the notch depth is about 10–15 mm (see Fig. 2 for size D= 300 mm). This is about 1.5 maximum aggregate size. The difference between the type 1 and type 2 fracture behaviors can also be seen in terms of the work of fracture. The work gddone in each strip (per unit ligament area) until the stress in that strip had been reduced to 0.75σt, was computed for each strip. However, many simulations (especially those with shallow notches) were terminated soon after the peak load was reached. In such cases, not many strips achieved the chosen stress limit of 0.75σt. The gdvalues from sufficiently damaged strips (averaged over the computed realizations, separately for each size and notch depth) are plotted in Fig. 3. The figure shows the crack initiation from a smooth surface to be followed by energy dissipation much larger than that for notched specimen. The energy comparison reveals no gradual transition. It is also interesting to compare the lattice-particle simulations to the standard cohesive 6 311 Jan Eli´aˇs and Zdenˇek P. Baˇzant Figure 4: Comparison of the relative nominal strengths obtained by the lattice-particle model, the cohesive crack model and the universal size effect law. crack modeling. The macroscopic cohesive law was obtained from stresses σxand energies gdmeasured in strips by the same procedure as introduced in [4]. For each size, a different macroscopic cohesive law was found by averaging the results from 6 realizations with α0=0.5. These cohesive laws were then approximated by piecewise linear functions. Then load–CMOD curves were computed for all the notch depths using the pseudo-boundary integral method [1] and piecewise linear approximations. The maximum loads obtained by the lattice-particle model are compared to the maximum loads from the cohesive crack model in Fig. 4 through relative nominal strength σN/σt. The nominal strength σNis defined as the maximal elastic stress in a notch-less beam loaded by peak force Fmax. Our beam geometry gives σN= 90Fmax/D. Assuming that the transition is a deviation from the cohesive crack model, one can see almost no transitional regime in this figure. Thus we see that the cohesive crack model is close to the lattice-particle model for most of the sizes and notch depths. Only very small sizes and very shallow notches lead to some deviations. However, as already pointed out, these notches have the depth of roughly the minimum aggregate radius and are thus poorly represented in the model. What is clear is that a discrepancy occurs at the zero notch depth, similar to what we concluded from the energy profiles in Fig. 3. 5 VARIATIONS IN DISSIPATED ENERGY The energy profiles in Fig. 3 show that the energy dissipated up to the point where the stress gets reduced to 0.75σtdepends on the specimen size. Whereas most simulations for D= 500 mm give gd≈20 N/m, the smaller specimens give lower values. Another systematic, though less visible, decrease of gdcan be seen for increasing notch depths. For the beam sizes of 500, 400 and 300 mm, the relative notch lengths α0=0.5 and 7 312 Jan Eli´aˇs and Zdenˇek P. Baˇzant Figure 5: Post-peak crack patterns when the load is decreased to 5% of the maximum load. 0.7 give a lower gdthan the shallower depths. This seems to be a continuous process originating from the width variation of the FPZ; i.e., the wider the zone, the more energy is dissipated. Visual comparisons of the damaged zone width can be seen in Fig. 5. The smaller the beam size and the deeper the notch, the narrower is the FPZ. Most of the energy is released in the middle of the zone and thus the energy profile comparison makes this dependency less obvious. The size dependence of gdis also projected into the reconstructed macroscopic cohesive law. The cohesive laws used in the pseudo-boundary integral method in the previous paragraph (and reconstructed for each size from the simulations with α0=0.5) are different. The larger the specimen size, the larger are the crack openings for a given cohesive stress. The conclusions from these numerical simulations broadly agree with experimental observations [10]. The fact that less energy is dissipated as the fracture process gets closer to the upper surface has already been pointed out [11]. Further it appears that the energy dissipation increases as the fracture process gets close to the bottom surface. 6 UNIVERSAL SIZE EFFECT LAW To describe the dependence of the nominal strength on both i) the relative notch depth α0and ii) the structure size, D, asymptotic matching is useful. The resulting formula, known as the universal size effect law, contained a discontinuity of slope [12, 13, 14, 15]. Here we consider an improved version from which the discontinuity has been removed [16, 17]; it reads: σN=EGf g 0cf+g0D1/21−rc2 fg 0e−kα2 0 4(lp+D)(g 0cf+g0D)1/r (1) where g0,g  0,g  0= values of the dimensionless LEFM energy release function evaluated at α=α0;g(α0) = square of dimensionless stress intensity factor and its derivative at α0;E 8 313 Jan Eli´aˇs and Zdenˇek P. Baˇzant Figure 6: Dependence of nominal strength σNon relative notch depth α0and specimen size Dobtained by a) lattice-particle simulations and b) analytical formula Eq. (1). = effective Young’s modulus; rand k= empirical positive parameters; lp= approximate width of the FPZ; cf= material constant, such that a=a0+cf= effective crack length at the peak load. Unlike the present model, which uses a deterministic material strength, the formula can also incorporate the effect of spatial variability of material strength which is the source of the statistical size effect. To determine the optimal parameters for the type 2 fracture, the Levenberg-Marquardt nonlinear optimization algorithm was used for least-square fitting of the average nominal strength values. Based only on the first part of Eq. (1) and on the nominal strength values for the relative notch depths larger than 0.05, the macroscopic initial fracture energy Gf= 30.3 N/m and cf= 19.3 mm was found. The value of Gfcorresponds to the initial fracture energy from the reconstructed cohesive law. The value of cfseems to be rather small, only 0.15lch, where lch =EGf/σ2 tis Irwin’s characteristic length. According to [1], cfshould be about π/24lch for rectangular softening and 0.419lch for linear softening. The remaining parameters lp,kand rwere found by including the shallower notches and the second part of Eq. (1). But a second least-square fitting yielded unrealistically small values of rand lp(both were virtually zero). As another approach, it was tried to use realistically chosen values lp=0.01 m, k= 300 and r=1/2. But neither approach was able to fit the nominal strengths provided by the lattice-particle model for zero notch depth. A better agreement for zero notch depth was obtained by considering the entire surface and optimizing all the five parameters (Gf,c f,l p,k,r) simultaneously. But such improvement of the fit for a zero notch depth leads to disagreement for the intermediate notch depths; i.e., incorrect Gfand cfis found. The two-parameter surface σN(D,α0), obtained by averaging the nominal strength values from simulations with the lattice-particle model, is shown in Fig. 6a, whereas the surface given by Eq. (1) for optimized parameter values is shown in Fig. 6b. The averaged simulated nominal strengths are reported in Tab. 1. The fit by Eq. (1) is also added to Fig. 4. Note that the σNvalues from the universal 9 314