Installation — business terrible  1 part
September 8th, 2015
The power of the microplane model is limited by the assumption of a kinematic (or static) constraint, which is a simplification of reality. To avoid this simplification, one needs to tackle the boundary value problem of the growth of many statistically unifonnly distributed cracks in an infinite elastic body. This problem is not as difficult as it seems. It appears possible that an approximate solution of this problem might once supersede the microplane model as the most realistic predictive approach to cracking damage.
The problem of calculating the macroscopic stiffness tensor of elastic materials intersected by various types of random or periodic crack systems has been systematically explored during the last two decades and effective methods such as the selfconsistent scheme (Budianski and O’Conell 1976; Iloenig 1978), the differential scheme (Roscoc 1952; Hashin 1988), or the MoriTanaka method (Mori and Tanaka 1973) have been developed. A serious limitation of these studies was that they dealt with cracks that neither propagate nor shorten (Fig. 14.5.1a). This means that, in the context of response of a material with growing damage illustrated by the curve in Fig. 14.5.1a, these formulations predict only the secant elastic moduli (such as Es in Fig. 14.5.1a). Such information does not suffice for calculating the response of a body with progressing damage due to cracking.
To calculate the response of a material with cracks that can grow or shorten, it is also necessary to determine the tangential moduli, exemplified by Et in Fig. 14.5.1b. Knowledge of such moduli makes it possible, for a given strain increment, to determine the inelastic stress drop 6ocr (Fig. 14.5.1b). This problem has recently been studied by Bazant and Prat (1995, 1997) and Prat and Bazant (1997). Its solution will now be briefly reviewed.
We consider a representative volume V of an elastic material containing on the microscale many cracks (microcracks). On the macroscale, we imagine the cracks to be smeared and the material to be represented by an approximately equivalent homogeneous continuum whose local deformation within the representative volume can be considered approximately homogeneous over the distance of several average crack sizes. Let є and a be the average strain tensor and average stress tensor within this representative volume. To obtain a simple formulation, we consider only circular cracks of effective radius a.
Consider the material to be intersected by N families of random cracks, labeled by subscripts /t — 1,2,..N. Each crack family may be characterized by its spatial orientation?7)(, its effective crack radius «,,, and the number of cracks in family ц per unit volume of the material. The compliance tensor may be considered as the function C — С(а;, аг,…, ajy; щ, гі2,… ,rijv). Approximate estimation of this function has been reviewed by Kachanov and coworkers (Kachanov 1992, 1993; Sayers and Kachanov 1991; Kachanov, Tsukrov and Shafiro 1994).
To obtain the tangential compliance tensor of the material on the macroscale, the cracks must be allowed to grow during the prescribed strain increment Ьє. This means that the energy release rate per unit length
Figure 14.5.1 Stressstrain curves and moduli: (a) effective secant moduli; (b) tangential moduli and inelastic stress decrement due to crack growth.
of the front edge of one crack must be equal to the given critical value R(atl) (or to the fracture energy Gy, in the case of LEFM). For the sake of simplicity, we will enforce the condition of criticality of cracks only in an overall (weak) sense, by assuming that the average overall energy release rate of all the cracks of one family within the representative volume equals their combined energy dissipation rate.
Our analysis will be restricted to the case when the number of cracks in each family is not allowed to change (бпц = 0), i. e., no new cracks are created and no existing cracks are allowed to close. This does not seem an overly restrictive assumption because a small enough crack has a negligible effect on the response and can always be assumed at the outset. Besides, concrete is full of microscopic cracks (or flaws) to begin with, and no macroscopic crack nucleates from a homogeneous material.
The incremental constitutive relation can be obtained by differentiation of Hooke’s law. It reads 6e — CSa + (dC/da^cr 6aM from which
U=1
where 6 denotes infinitesimal variations and C(a(I) is the fourthorder macroscopic secant compliance tensor of the material with the cracks, and E(aJt) is the fourthorder secant stiffness tensor, whose 6×6 matrix is inverse of the matrix of C(a;i).
The surface area of one circular crack of radius o. lL is A — and the area change when the crack radius increases by 6a is 6AM — 2ттаІІ6а1і. We assume we can replace the actual crack radii by their effective radius aд.
The crack radius increments 6aд must be determined in conformity to the laws of fracture mechanics. Let us assume that the cracks (actually microcracks) can be described by equivalent linear elastic fracture mechanics (LEFM) with an Rcurve R(atl). This means that the energy release rates must be equal to (or Gf, in the case of LEFM). For the special case Л(ам) — Gj = fracture energy of the matrix of the material, the cracks follow LEFM.
The complementary energy density of the material is U = U (a, aM) = cr ■ С(aM)cr, where the dot indicates scalar product of two secondorder tensors. To make the problem tractable, we impose the energy criterion of fracture mechanics (energy balance condition) only in the overall, weak sense. This leads to the following N conditions of crack growth (Bazant and Prat 1995; Prat and Bazant 1997):
(repetition of subscript/tin this and subsequent equations does not imply summation); — crack growth indicator which is equal to 1 if the. crack is growing (6atl > 0) and 0 if it is closing (Safl < 0), while any value 0 < kil < 1 can be used if 6aд ~ 0.
Differentiation of (14.5.2) provides the incremental energy balance conditions;
dC – A, ( QLC
cr. 6a + ( – jer • Де Qa~a " ^п^КцСуб^и J 6a„ — 2m^atlGy6(/t — 1,2,…N)
(14.5.3)
where бкц — 0 except when the crack growth changes to no growth or to closing, or vice versa. The handling of the large jump in ajL is exact if — s” and atl — a“ld because 6(K, tlatl) — – (ісдам)°и – – І exactly.
Substitution of (14.5.1) into (14.5.3) leads to a system of N equations for N unknowns 6ci….6ajy:
N
J2A = (F=l,N) (14.5.4)
u=
where
В/г ~ & ■ 7Г—E6e — 27rn/ta°ldG/бкр (14.5.6)
0(1^
After solving (14.5.4), evaluating 6cr from (14.5.1) and incrementing cr, one must check whether the case 6afl > 0 and ONp < 0 (or ejy < 0) occurs for any /t. If it does, the corresponding equation in the system (14.5.4) must be replaced by the equation 6ap = 0. The solution of such a modified equation system must be iterated until the case 6atl > 0 and < 0 (or cfy < 0) would no longer occur for any /л.
The formulation needs to be further supplemented by a check for compression. The reason is that the energy expression is quadratic, which means that (14.5.2) is invariant when cr is replaced by —cr. Thus, a negative stress intensity factor K/lL can occur even when (14.5.4) is satisfied, and so the sign of K/p must be checked for each crack family )i. The case K/p < 0 is inadmissible.
Since the present formulation yields only the values of (Kip)2 ~ ETZ(at,) but not the values of Kip, the sign of Kip must be inferred approximately. It can be considered the same as the sign of the stress ONp — iip ■ сгтір in the direction normal to the cracks of 71th family («,, = unit vector normal to the cracks). [Alternatively, the sign of Kip could be inferred from the sign of the normal component eiv> = "a ‘ £л ”п cracking strain tensor ecr = CCTcr. J Approximate though such estimation
surely is on the microscale of an individual crack, it nevertheless is fully consistent with the macroscopic approximation of C implied in this model. The reason is that all the composite material models for cracked solids are based on the solution of one crack in an infinite solid, for which the sign of Kip coincides with the sign of ацр (or ).
Usually the six independent components of Se are known or assumed, and then (14.5.4) represents a separate system of only N equations for Sai, …бац (this is a. simplification compared to the formulation in the paper, which led to a system of N 1 6 equations). In each iteration of each loading step, the values of Kp are set according to the sign of 6atl in the preceding step or preceding iteration.
If Sap — 0, or ifi(due to numerical error) ба^ is nonzero but less than a certain chosen very small positive number 6, Up is arbitrary and can be anywhere between 0 and 1, which makes equation (14.5.3) superfluous. The condition (14.5.2) of energy balance in the of constant crack length becomes meaningless and must be dropped. It needs to be replaced by an equation giving Up (or Kip) as a function of cr), which must be used to check whether Kp indeed remains within the interval (0,1). However, it seems that for proportionally increasing loads the case Sap — 0 is not important and its programming could be skipped, using conveniently the value Kp ~ 0.5.
The foregoing solution was outlined in Bazant and Prat (1995) and worked out in detail in Prat and Bazant (1997) (with Addendum in a later issue).
To be able to use (14.5.5), we must have the means to evaluate the effective secant stiffness C as a function of ttp. Bazant and Prat (1995) and Prat and Bazant (1997) adopted the approximate approach developed by Sayers and Kachanov (1991) using the symmetric secondorder crack density tensor
N
OL = nval ® % (14.5.7)
д=1
(Vakulenko and Kachanov 1971; Kachanov 1980, 1987b). In this approach, the effective secant compliance C is derived from an elastic potential F which is considered as a function of (he crack density tensor cr (in addition of the stress tensor cr):
The elastic potential F(a, cr) can be expanded into atensorial power series. Sayers and Kachanov (1991) proposed to approximate potential F by a tensor polynomial that is quadratic in a and linear in cr:
in which 1)1 and i]2 are assumed to depend only on a — tr a = ^2 rnap’ l’,c *’rsl invariant of a (Sayers and Kachanov 1991). The strain tensor follows from (14.5.9):
The functions rji (a) and 772(a) can be obtained by taking the particular form of the preceding formulation for the case of random isotropically distributed cracks and equating the results to those obtained using, e. g., the differential scheme (Ilashin 1988). To this end, we note that if the orientation, density, and size of cracks is isotropically distributed, then a must be spherical, and since its trace is equal to a, it must be a = (a/3)l. Substituting this into the preceding equation we get for this case:
Noting now that the resulting behavior in (14.5.11) must be elastic with effective elastic modulus Ecff and effective Poisson’s ratio Ум, we get the following relationship between the functions 77;, the effective elastic constants and the a:
where E and у are the elastic constants of the uncracked material (included in C°); here we indicate explicitly that the effective elastic constants depend on ft. This dependence can be obtained, for example, by using the differential scheme which, for quasibritlle materials, gives better predictions than the self – consistent scheme (Bazant and Prat 1997). The resulting relationships (see e. g., Hashin 1988 for the details of the derivation) are as follows:
From Eqs. (14.5.7), (14.5.10), and (14.5.12)—(14.5.14) the effective compliance tensor is obtained as a function of aM. Then Saй is calculated from F. q. 14.5.4 for a given be as indicated before.
The crack density may be characterized as a continuous function nlL of spherical angles ф and 0 (Prat and Bazant 1997). Function ntl is then sampled at spherical angles ф^ and (7fI corresponding to the orientations of the microplanes in the microplane model. For isotropic materials such as concrete, the distribution of is initially uniform, and a very small but nonzero value must be assigned to every as the initial condition because no new cracks are allowed to nucleate. For an initially anisotropic material such a rock with joints, function is assumed to peak at a few specified discrete orientations ф*,в*.
In ongoing studies, the //curves are used by Bazant and coworkers as a means to approximate the effect of plastic strains in the matrix of the material occurring simultaneously with the crack growth. Unlike classical plasticity, the plastic strain cannot be considered here to be smeared in a continuum manner because the cracks cause stress concentrations. Therefore, plasticity of the matrix will get manifested by the formation of a plastic zone at the front edge of each microcrack. As is well known, the effect of such a plastic zone can be approximately described by an Яcurve.
The lattice models or particle systems are computationally very demanding and a very efficient computational algorithm must be used. A highly efficient algorithm, which was applied to simulation of sea ice fracture, is presented in Jirasek and Bazant (1995a). It is an explicit algorithm for fracture dynamics, but it can also be used for static analysis in the sense of dynamic relaxation method. (Fracturing in particle systems with more than 120,000 degreesoffreedom was simulated with this algorithm on a desktop 1992 work station.) Computational effectiveness will be particularly important for three dimensional lattices, the use of which is inevitable for obtaining a fully realistic, predictive model.
As it now stands, the lattice or particle models can provide a realistic picture of tensile cracking in concrete in twodimensional situations. However, solution of significant threedimensional problems or nonlinear triaxial behavior as well as simulation of behavior in which compression and shear fracturing is important is still beyond reach. Thus, the lattice or particle models have not attained the degree of generality already available with the finite element approach.
Although computer programs for lattices are attractive by their simplicity, it must nevertheless be recognized that a lattice modeling of a continuum is far less powerful than the finite element method because the stress and strain tensors cannot be simulated by the elements of the lattice and, thus, nonlinear tensorial behavior cannot be directly described. From this viewpoint, the numerical concrete of Roelfstra and Wittmann, seems preferable to the lattice model of van Mier and Schlangen because the particles and mortar are discretized by a much finer mesh of finite elements.
To sum up, lattice models or particle systems have proven to be a useful tool for understanding fracture process’and clarifying some relationship between the micro – and macro – characteristics of quasibrittle heterogeneous materials. These models appear to be particularly suited for failures due to tensile fracturing
and capture well the distributed nature of such fracturing and its localization. However, one must keep in mind that these models, in their present form, cannot simulate threedimensional situations, larger structures (even twodimensional), compressive and shear fractures, and nonlinear triaxial stressstrain relations. Overall, these models are still far inferior to finite element models and do not really have a predictive capability. No doubt significant improvements may be expected in the near future.
Bazant, Tabbara et al. (1990) used the random particle system described before (Figs. 14.4.lab and 14.4.2) to simulate tensile tests and bending tests on notched specimens. A similar model was used by Jirasek and Bazant (1995a, b) to relate the microscopic features of the model (such as the softening curve and the statistics of strength distribution) to the macroscopic properties, particularly size effect and fracture energy.
Fig. 14.4.6 shows direct tension specimens of various sizes studied computationally by Bazant, Tabbara et al. (1990), with the results displayed in Fig. 14.4.7 as the calculated curves ofload (axial force resultant)
Figure 14.4.6 Geometrically similar specimens of various sizes with randomly generated particles (adapted from Bazant, Tabbara et al. 1990). 
vs. relative displacement between the ends.
Fig. 14.4.7b gives the curves for several specimens of the smallest size from Fig. 14.4.6. Fig. 14.4.7c shows the curves for the medium size specimens and Fig. 14.4.7d the curve for one large size specimen. Fig. 14.4.7e shows, in relative coordinates, the average response curves calculated for the small, medium, and large specimens. Note that while the prepeak shape of the load displacement curve is size independent, the postpeak response curve is getting steeper with increasing size.
Fig. 14.4.8 shows the progressive spread ol’cracking in one of the smallest specimens from Fig. 14.4.6. The cracking patterns arc shown lor four different points on the load displacement diagram, as seen in Fig. 14.4.8a, the first point corresponding to the peak load. The dashed black lines are the normals to the links that undergo softening and correspond to partially formed cracks. The solid lines are normal to completely broken links and represent fully formed cracks. The gray dashed lines represent normals to the links that partly softened and then unloaded, and correspond to partially formed cracks that are closing. Note from Fig. 14.4.8 that the cracking is at first widely distributed, but then it progressively localizes.
Fig. 14.4.9 shows the calculated peak loads for the specimens from Fig. 14.4.6 in the usual size effect plot of the logarithm of nominal strength vs. logarithm of the size. Bazant, Tabbara et al. (1990) interpreted the results in terms of the classical siz. e effect law Eq. (1.4.10) with relatively good results. The recent results of Bazant explained in Section 9.1, particularly the size effect formula for failures at crack initiation from a smooth surface (Section 9.1.6) suggest that these results must be interpreted using Eq. 9.1.42. Thus the results of Bazant, Tabbara et al. (1990) have been fitted here by the simplest version of this curve (for which 7 — 0). Fig. 14.4.9 shows the resulting fit, which is excellent for the mean, values of the data. .
Threepointbend fracture specimens of three sizes in the ratio 1:2:4 wete simulated in the manner illustrated in Fig. 14.4.10a. Fig. 14.4.10b shows the size effect plot obtained from the three sizes of the specimens in Fig. 14.4.10a for three different materials. As can be seen, the calculated maximum loads can be well approximated by Bazant’s size effect law, Eq. (1.4.10).
By fitting the size effect law to the maximum load obtained by the lattice or particle model for similar specimens of different sizes, one can determine the macroscopic fracture energy C j of the particle system and the effective length Cf of the fracture process zone (see Chapter 6). It thus appears that fracture simulations with the lattice model or random particle system provide a further verification of the general applicability of the size effect law.







Figure 14.4.10 Simulation of three threepointbend tests by random particle model.(a) Thrccpointbend specimens with d = 36.72, and 144 mm and (b) corresponding size effect plot. (Adapted from BaiSant, Tabbara et al. 1990.)
At the same time, the size effect law is seen to be an effective approach for studying the relationship between the microscopic characteristics of the particle system, simulating the microscopic properties of the material, and the macroscopic fracture characteristics.
Such studies have been undertaken by Jirasek and Bazant (1995b). Fig. 14.4.11 show the results of a large number of such simulations, dealing with two dimensional threepointbend fracture specimens of different sizes. In these specimens, the microductility number, representing the ratio sc/ep, was varied (see Fig. 14.4.2b). The coefficient of variation of the microstrength of the particle length, used in random generation of the properties of the links, was also varied (the microstrength was assumed to have a normal distribution).
It was found that both the microductility and the coefficient of the microstrength of the links have a significant effect 011 the macroscopic fracture energy Gy and on the effective length (.7 of the process zone; see Fig. 14.4.1 la, c. Randomness of these plots is largely due to the fact that the number of simulations was not very large (the values in the plot are the averages of the values obtained in individual sets of simulations of specimens of different sizes). The plots in Fig. 14.4.1 la, c have been smoothed as shown in Fig. 14.4.1 lb, d by the following bilinear polynomials which provide optimum (its:
Figure 14.4.11 Normalized fracture energy and normalized effective process zone size as a function of two parameters: (a) and (c) computed, (b) and (d) fitted by bilinear functions (from Jirdsek and Baz. ant 1995b). 
Cf = 0.64 – I – 0.08wy + 0.097/ — 0.19u>/7/, (14.4.2)
in which the superimposed bars refer to average values, and cjj and 7/ are the coefficient of variation of microstrength and the microductility number. Obviously, the effect of various other microscopic characteristics of particle systems on their macroscopic fracture properties could also be studied in this manner, exploiting the size effect law.
An important aspect of the model is the generation of the lattice configuration. In many works regular lattices have been used. However, recently Jirasek and Bazant (1995a), and also Schlangen (1995)
Figure 14.4.5 Failure patterns for various values of a: (a) 0°, (b) 22.5°, (c) 45° using a regular lattice with random strength, elastic stiffness and microl’racture energy of the links (from Jirasek and Bazant 1995a).
demonstrated that a regular lattice always impresses a strong bias on the direction of fracture propagation.
For the square lattices with diagonals analyzed by Jirasek and Bazant, it is, of course, possible to choose the elastic stiffnesses of the links in the main directions of the square mesh and the diagonal directions, the corresponding strength limits of the links and the corresponding microfraclure energies in such ratios that the lattice is isotropic in terms of elastic properties, strength along straight line cuts, and fracture energies dissipated on such cuts for any orientation of the cut. However, even in that case, the fracture tends to run preferentially among the mesh lines. This has been blatantly demonstrated by simulations of fracture of a circular specimen on which a regular square mesh with diagonals was overlaid; see Fig. 14.4.4. In this particle simulation the fracture was caused by an impact at the bottom of the circle in upward direction. In Fig. 14.4.4a the impact was in the direction of the square mesh lines, in Fig. 14.4.4d in the direction of the diagonals, and in Fig. 14.4.4b and 14.4.4c in two intermediate directions. Note the enormous differences in fracture patterns, which were also manifested by great differences in peak loads and energies dissipated. When all the properties of the links of a regular lattice were randomized, strong directional bias of fracture still remained; see Fig. 14.4.5.
Only when a geometrically random lattice was used in Jirasek and Bazanl’s (1995a) study, the directional bias was eliminated, except for small random differences between meshes. Similar results were found by Schlangen (1995) for a doubleedge notched specimen subjected to shear. These results indicate that random (unstructured) lattices must be used to avoid directional bias.
The simplest model is a pinjointed truss, in which only the centertocenter forces between the particles are considered (Fig. 14.4.lab, Bazant, Tabbara et al. 1990). A more refined model is that ofZubelewicz. and Bazant (1987), which imagines rigid particles separated by deformable thin contact layers of matrix that respond primarily by thickness extensioncontraction and shear (Fig. 14.4.1 cd). Since the intemodal links also transmit shear, moment equilibrium of the nodes needs to be considered, while for the pinjointed truss it need not. Therefore, this model has three degreesoffreedom per node (two displacements and one rotation, with corresponding two force components and one moment) for planar lattices, while the pinjointed model has only two degre. esoffreedom per node. In the spatial case, the model ofZubelewicz and Bazant requires six degreesoffreedom per node, i. e., three displacements and three rotations, while the pinjointed The simplest model truss requires only three degreesoffreedom per node. There is an additional important advantage of shear transmissionThe simplest model —it makes it possible to obtain with the lattice any Poisson ratio, while a random or regular pinjointed lattice (truss) has Poisson ratio always 1/3 in two dimensions and 1/4 in three dimensions.
The simplest model In the model of Bazant, Tabbara et al. (1990) and Zubelewicz. and Bazant (1987), the major particles in the material (large aggregate pieces) are imagined as circular and interacting through links as shown in Fig. 14.4.2. In the initial work ofZubelewicz. and Bazant, the link between particles was assumed to transmit both axial forces and shear forces, the latter based on the rotations of particles. In the subsequent model by Bazant, Tabbara et al. (1990), the particle rotations and transmission of shear were neglected and only axial forces were assumed to be transmitted through the links. In such a case, the system of particle links is equivalent to a truss. As pointed out before, the penalty to pay for this simplification is that the Poisson ratio of a random planar truss is always 1/3 (and for a spatial truss 1/4). Another consequence of ignoring particle rotations and interparticle shears is that the fracture process zone obtained becomes narrower. But this can be counteracted by assuming a smaller postpeak softening slope for the interparticle stress displacement law, and also by introducing a greater random scatter in the link properties, both of which tend to widen the fracture process zone.
A random particle configuration must be statistically homogeneous and isotropic on the macroscale. In the simulation of concrete, the configuration must meet the required granulometric distribution of the particles of various sizes, as prescribed for the mix of concrete. The problem of generation of random configurations of particles in contact under such constrains involves some difficult and sophisticated aspects (see, e. g., Plesha and Aifantis 1983).
However, the problem becomes much simpler when the particles do not have to be in contact, as is the case for aggregate pieces in concrete. In that case, a rather simple procedure (Bazant, Tabbara et al. 1990) can proceed as follows: (1) using a random number generator, coordinate pairs of particle centers (nodes)
arc generated one after another, assuming a uniform probability distribution of the coordinates within the area of the specimens; (2) for each generated pair a check for possible overlaps of the particles is made, and if the generated particle overlaps with some previously generated one, it is rejected; (3) the random generation of coordinate pairs proceeds until the last particle of the largest size has been placed within the specimen, (4) then the entire random placement process is repeated for the particles of the next smallest size, and then again for the next smallest size, etc. (The number of particles of each size is determined in advance according to the prescribed mix ratio and granulometry.)
To determine which particles interact, a circle of radius i>r, is drawn around each particle i (with ft ~ 5/3) as shown by the dashed lines in Fig. 14.4.2a: Two particles interact if their dashed circles intersect each other. See Bazant, Tabbara et al. (1990) for the details of the assignment of the dimensions of the truss element, particularly the crosssection Am and length Lm of the deformable portion (labeled with subscript m for matrix). In a later study by Jirasek and Bazant (1995a), a uniform stiffness of all the links was assumed.
Fig. 14.4.2b shows a typical computergenerated random particle arrangement resembling concrete, and the corresponding truss (random lattice). Fig. 14.4.2c shows the stressstrain relation for the interparlicle links, characterized by the elastic modulus Em, tensile strength limit /4m, and the postpeak softening slope Es (or alternatively by strain су at complete failure, or by G™, each of which is related to the foregoing three parameters). The microscopic fracture energy of the material, G™, is represented by the area under the stressstrain curve in Fig. 14.4.2c, multiplied by the length of the link. The ratio of є/ to the strain єр at the peak stress may be regarded as the microductility of the material.
The lattices in Fig. 14.4.1a–d attempt to directly simulate the major inhomogeneities in the microstructure of concrete. By contrast, the model introduced by Schlangen and van Mier (1992) takes a lattice (in the early versions regular, but later randomized) that is much finer than the major inhomogeneities. Its nodal locations and links are not really reflections of the actual microstructure (Fig. 14.4. lef). Rather, the microstructure is simulated by giving various links different properties, which is done according to the match of the lattice to a picture of a typical aggregate arrangement.
Van Mier and Schlangen take advantage of the available simple computer programs for frames and assume the lattice to consist of beams which resist not only axial forces but also bending. Due to bending, the internodal links (beams), of course, also transmit shear, same as in the model of Zubelewicz and Bazant (Figs. 14.4.1cd and ef). This feature is useful, because shears are indeed transmitted between adjacent aggregate pieces and across weak interfaces in concrete, and because arbitrary control of the Poisson ratio is possible. However, the idea of bending of beams is a farfetched idealization that has nothing to do with reality. No clear instances of bending in the microstructures of concrete can be identified.
The idealization of the links as beams subject to bending implies that a bending moment applied at one node is transmitted to the adjacent node with the сапуover factor 0.5, as is well known from the theory of frames. This value of the carryover factor is arbitrary and cannot not have anything in common with real behavior. In the model of Zubelewicz. and Bazant (Fig. 14.4.1 cd), the shear resistance also causes a transmission of moments from node to node, however, the carryover factor is not 0.5 and can have different values. The transmission of moments is, in that model, due to shear in contact layers between particles, which is a clearly identifiable mechanism. In consequence of this analysis, it would seem better to consider the carryover factor in the lattices of van. Mier and Schlangen to be an arbitrary number, determined either empirically or by some microstructural analysis. This means that the 6×6 stiffness matrix for the element of the lattice, relating the 6 generalized displacement and force components of a beam sketched in Fig. 14.4.If, should be considered to have general values in its off diagonal members, not based on the bending solutions for a beam but on other considerations. In fact, the use of such a stiffness matrix would require only an elementary change in the computer program for a frame (or lattice with bending). Of course, if the need for such a modification is recognized, the model and van Mier and Schlangen becomes essentially equivalent to that of Zubelewicz and Bazant, except that the nodes do not represent actual particles and the lattice is much finer than the particles.
The beams’in the lattice model of Schlangen and van Mier are assumed to be elasticbrittle, and so, when the failure criterion is met at one of the beams, the link may be removed. This means that at each step the computation is purely elastic. This is computationally efficient, but makes the model predict a far too brittle behavior, even for threedimensional lattices (van Mier, Vervuurt and Schlangen 1994).
Another important aspect in lattice models is the size of the links. Unlike the lattice of Zubelewicz and Bazant, which directly reflects the particle configurations and thus cannot (and should not) be refined, the
Figure 14.4.3 Dependence of the loaddisplacement curve on the size of Ihe lattice links: (a) basic element geometry, (b) loaddisplacement curves for various lattice spacings (adapted from Schlangen 1995).
lattice of van Mier and Schlangen lias an undetermined nodal spacing, which raises additional questions. First, as is well known, frames or lattices with bending are on a large scale asymptotically equivalent to the socalled micropolar (or Cosserat) continuum (e. g., Baz. ant and Ccdolin 1991, Sec. 2.102.11). Pinjointed trusses, on the other hand, asymptotically approach a regular continuum on a large scale. The micropolar continuum is a continuum with nonsymmetric shear stresses and with couple stresses. It possesses a characteristic length, which is essentially proportional to the typical nodal spacing of the lattice approximated by the micropolar continuum. While, in principle, the presence of a characteristic length is a correct property for a model of concrete, the characteristic length should not be arbitrary but should be of the order of the spacing of the major aggregate pieces. In this regard, the lattice of van Mier and Schlangen appears to be too refined. Moreover, as transpired from recent researches and the previous chapter on nonlocal concepts, the micropolar character or the presence of characteristic length should refer only to the fracturing behavior and not to the elastic part of its bonds. The model of Schlangen and van Mier goes against this conclusion, since even the elastic response of the lattice is asymptotically approximated by a micropolar rather than regular continuum.
Furthermore, a question arises about the dependence of the response on the lattice spacing. A recent study of Schlangen (1995) shows that the crack pattern is not strongly affected by the size of the beams, but the loaddisplacement is affected in much the same way as mesh refinement in local strainsoftening models: the finer the lattice, the less the inelastic displacement and the dissipated energy, as illustrated in the loadcrack opening curves in Fig. 14.4.3 for a square specimen subjected to pure tension. Indeed, it is easy to imagine that upon infinite refinement the stresses in a beam close to a crack lip must tend to infinity and thus a precracked specimen must fail for a vanishingly small load (roughly proportional to the square root of the beam siz. e). Also, the shorter the beams, the smaller the dissipated energy, because the volume of material affected by the crack is smaller the smaller the elements. Note that the lattice analyzed by Schlangen in Fig. 14.4.3 has random strength in all the cases with identical probabilistic distribution. Therefore, randomness does not relieve mesh sensitivity as sometimes claimed.
The meshsensitivity of Schlangen and van Mier’s model can probably be artificially alleviated or eliminated by taking a beam strength inversely proportional to the square root of the beam size, similar to the equivalent strength method described in Section 8.6.4 for crack band analysis. However, this is purely speculative and a more sound basis should be built for the lattice models before they can be confidently used as predictive (rather than just descriptive) models. A nonlocal fracture criterion may serve as an alternative solution to the problem, but this would break the computational efficiency of the elasticbrittle beam lattice model. Note that nonlocality (i. e., interaction at finite distance) is automatically implemented in the particle models because the particle distances are finite and fixed, so the lattice siz. e should also be fixed as dictated by the microstruclurc.
A large amount of research, propitiated by the advent of powerful computers, has been devoted to the simulation of material behavior based directly on a realistic but simplified modeling of the microstructure —its particles, phases, and the bonds between them. A spectrum of diverse approaches can be found in the literature spanning an almost continuous transition from the finite element simulations, with the classical hypothesis of continuum mechanics, to discrete particle models and lattice models in which the continuum is approximated a priori by a system of discrete elements: particles, trusses, or frames.
An extreme example of the continuum approach —in view of the fineness of material subdivision— is the numerical concrete of Roclfstra, Sadouki and Wittmann (1985), Witlmann, Roelfstra and Kamp (1988) and Roelfstra (1988), in which the mortar, the aggregates, and their interfaces are independently modeled by finite elements. This requires generation of the geometry of the material (random placement of aggregates within the mortar) and the detailed discretization of the elements to adequately reproduce the geometry of the interfaces. With a completely different purpose, but with the same kind of analysis, Rossi and Richer (1987) and Rossi and Wu (1992) developed a random finite element model in which the microstruclure is not directly modeled, but is taken into account by assigning random properties to the element interfaces. The common feature of these approaches is that, before cracking starts, the displacement field is approximated by a continuous function.
The particle and lattice models do not model the material continuously, but substitute the continuum by an array of discrete elements in the form of particles in contact, trusses, or frames, in such a manner that the displacements are defined only at the centers of the particles, or at the nodes of the truss or frame.
The origin of the particle approach can generally be traced to the development of the socalled distinct element method by Cundall (1971, 1978), Serrano and RodriguezOrtiz (1973), RodriguczOrtiz (1974), Kawai (1980), and Cundall and Struck (1979) in which the behavior of particulate materials (originally just cohesionless soils and rock blocks) was analyzed simulating the interactions of the particles in contact. This kind of analysis, which deals with a genuine problem of discrete particle systems, used highly simplified contact interaction laws permitted by the fact that the overall response is controlled mainly by kinematic restrictions (grain interlock) rather than by the details of the forcedeformation relation at the contacts. However, although the kinematics of the simulations appeared very realistic, the quantitative stressstrain (averaged) response was not quite close to the actual behavior. This shortcoming, which still
persists in many modem particle and lattice models, is largely caused by the fact that the simulations are usually twodimensional while a realistic simulation ought to be threedimensional.
The basic idea of the particle model can be extended to simulate the particular structure of composite materials, for example, the configuration of the large aggregate pieces of concrete, as done by Zubelewicz (1980,1983), Zubelewicz and Mroz(1983), and Zubelewicz and BaXant (1987), or the grains in a rock (Ple – sha and Aifantis 1983). In these cases the model requires defining the force interaction between particles (aggregates or grains) which are caused mainly by the relative displacements and rotations of neighboring particles. Although, for computational purposes, the problem is reduced to a truss (Fig. 14.4. lab) or to a frame (Fig. 14.4. lcd), the basic ingredient of such models is that the geometry (size) of the truss or frame elements and their properties (stiffness, strength, etc.) arc dictated by the geometry of the physical structure of the material (stiffness, size, shape, and relative position of aggregates or grains).
In contrast to this, the pure lattice models replace the actual material by a truss or frame whose geometry and element sizes are not related to the actual internal geometry of the material, but are selected freely by the analyst. The truss approach to elasticity, elementary atomistic representations of the physics of elasticity (i. e., arrays of atoms linked by springs shown in textbooks of solid state physics), was already proposed as early as 1941 by Hrennikoff. The lattice models have been championed by theoretical physicists for the simulation of fracturing in disordered materials (Herrmann, Hansen and Roux 1989; Channel, Roux and Guyon 1990; Herrmann and Roux 1990; Herrmann 1991) and have been developed to analyze concrete fracture by Schlangen and van Mier at Delft University of Technology (Schlangen and van Mier 1992; Schlangen 1993, 1995; van Mier, Vervuurt and Schlangen 1994). In their approach, a regular triangular frame of side length less than the dimensions of the smallest aggregates, is laid over the actual material structure (Fig. 14.4.1 e—F) and the properties of each beam are assigned according to the material the beam lies over, mortar, aggregate or interface. However, to eliminate directional bias of fracture, the lattice must be random (see Section 14.4.2).
This section presents a brief overview of the main concepts and results of the particle and lattice models as far as concrete fracture is concerned.
Figure 14.4.2 Random particle model of Bazant, Jirasek et al. (1994): (a) two adjacent circular particles with radii ті and and corresponding truss member ij; (b) typical randomly generated specimen and its corresponding mesh of truss elements; (c) constitutive law for matrix. (Adapted from Bazant, Tabbara et al. 1990.)
In unconfined straining, the microplane model displays softening. Therefore, localization limiters of some kind must be used to avoid spurious localization and nresh sensitivity, as for all other models with strain softening. This can be easily implemented using a nonlocal adaptation of the microplane model in which the inelastic stress increment is made nonlocal following the theory of microcrack interactions presented in the previous chapter (Bazant 1994b; see §13.3). This approach affects the flow of calculation only partially and a general finite element scheme can be used. Fig. 14.1.1c shows the basic calculation flowchart for this approach in which the nonlocal adaptation is implemented just after the microplane stresses get computed; the flow bifurcates and the inelastic incremental stress is computed following the nonlocal theory with microcrack interactions.
The microplane model as presented, or for that matter any constitutive model for damage, gives a prescription to calculate the stress tensor or as some tensorvalued function It of the strain tensor є (and of some further parameters depending on the loading history, e. g., on whether there is loading or unloading). So, o’ = Іі(є). The most robust (although not always the most accurate) method of structural analysis is to base the solution of a loading step or lime step on the incremental elastic stressstrain relation with inelastic strain involving the initial elastic moduli tensor E, as explained in Section 13.3.1. Then, for a local formulation, the inelastic stress increment tensor AS defined in (13.3.1) can be computed as
AS – E( ^IlCW £old) — ІЦЄпет/) + ЩЄоІсі) (14.3.1)
in which subscripts old and new label the old and new value of the variables at the beginning and end of the loading step (or time step); and S is the inelastic stress tensor due to nonlinear behavior. This stressstrain relation is used for both dynamic explicit analysis and static implicit analysis (as the iterative initial stiffness method).
A possible simple approach to introduce nonlocal effects is similar to the isotropic scalar nonlocal approach (PijaudierCabot and Bazant 1987), which was applied to the microplane model by Bazant and Ozbolt (1990, 1992) andOz. bolt and Bazant (1992). In this approach, the elastic parts of stress increments are calculated locally. The inelastic parts of the increments of S must be calculated nonlocally. This is accomplished by first determining, at each integration point of each finite, element, the average (or nonlocal) strains є, and then calculating nonlocal AS from these, i. e.,
Ac — ЕАє – AS; AS = E(WM – £„c„) – R(eIK. w) + К(є0м) (14.3.2)
The only modification required in a local finite element program is to insert the spatial averaging subroutine just before the calculation of AS.
A better approach is to introduce the crack interaction concept explained in Section 13.3 and write the incremental elastic stressstrain relation as
Дог =■ ЕАє – AS (14.3.3)
in which AS is given by Eq. (13.3.9). The spatially averaged strains are not calculated in this approach. The nonlocal part of the analysis proceeds in the following steps:
1. First, AS is calculated (in the local form) from (14.3.1) according to the microplane model. Then one calculates at each integration point of each finite element the maximum principal direction vectors n(4 [і = 1,2, 3) of strain tensor e, for which the value of e0id may be used as an approximation.
Figure 14.3.1 Local representative volume: orientation and size (adapted from Ozbolt and Bazant 1996). 
2. Then one starts a loop on principal strain directions n‘ (і — 1,2,3) of tensor є and evaluates
the inelastic stress changes in the directions that is, ASM = ■ ASn^ or A=
ASijn, j’ For those principal directions for which Дб’М < 0, the nonlocal calculations are skipped because the inelastic strain is not due to cracking, i. e., one jumps directly to the end of this loop; here we make the assumption that, on the microscale (but not on the macroscale), there is no softening in compression, which is true for the microplane model.
3. The values of ASW for the integration points of finite elements are then spatially averaged:
A“4l) = ~ it Д5<4Ч„Д V„ (14.3.4)
V It,
‘ V— 1
where V’fj — Yl*=i AV„ = normalizing factor, n = number of all the integration points inside the averaging volume, and afiu = given weight coefficients, whose distribution is suitably chosen with a bell shape in both x and у directions, described by a polynomial of the fourth degree. The bell shape, which is similar to that in the nonlocal damage approach (Chapter 13, Eq. (13.1.5)) is reasonable in that it gives larger contributions to the sum from points that lie closer. Because the spacing of major cracks in concrete is approximately the same as the spacing of the largest aggregate pieces, the size of the averaging volume may be assumed to be approximately proportional to the maximum aggregate size, da. For twodimensional analysis, the region of averaging should probably be taken as a rectangle with its longer side in the direction normal to (Fig. 14.3.1)
4. The values of the nonlocal principal inelastic stress increments Amust then be solved from the system of linear equations (13.3.31) based on the crack influence function However, as discussed before, exact solution is normally not needed. Depending on the type of program, one of two approximate methods can be used:
(a) In programs in which the loading step is iterated, these equations may be solved iteratively within the same iteration loop as that used to solve the nonlinear constitutive relation, using the following equation:
ASfnew ASW I — ДКЛ^Л5»оШ (м – 1,2, …N) (14.3.5)
in which A’ = crack influence matrix defined in Eq. (13.3.43), which must, however, be adjusted with factor kb for finite elements close to the boundary of concrete (Section 13.3.9).
(b) In explicit finite element programs without iteration, one may calculate from (14.3.5) only the first iterate (r — 1), which represents one explicit calculation, requiring only the values of Aand (Д5,?). The premise of this approximation is that the repetitions of similar calculations (for r — 1) in the next loading step (or time step) effectively serve as the subsequent iterations (for r = 2, 3,4,…) because the loading steps in the explicit programs are very small. This of course means that the correct value of А6’,^ gets established with a delay of several steps or time intervals (in other words, the computer program is using nonlocal inelastic stress increments that are several steps old; the nonlocal interactions expressed by the crack influence function are delayed by several steps).
We recall from Section 13.3.9 that the adjustment by factor кь must ensure that, even if part of the influencing volume protrudes beyond the boundary, the condition А’/ш 0 be met. Because
this condition may be written as ^interior ^ /^boundary ~ 0>lllc following adjustment is needed for the integration points of the elements adjoining the boundary;
= ki, A,lL,, кь — – Y, AfW / Y, Л,,„ (14.3.6)
interior a ‘ boundary ;/
For the remaining integration points in the interior, no adjustment is done, i. e., А!,ш — Л.
5. At each integration point of each finite element, the nonlocal inelastic stress increment tensor is then constituted from its principal values according to the following equation:
з з
ASkl = ]T AS^n^n^ or AS = Y 5(!>n(,) Jo n(0 (14.3.7)
І~. 1 І–1
based on the spectral decomposition theorem of a tensor.
Note that if, at some integration point, all the principal values of tensor AS arc nonpositive, then the foregoing nonlocal procedure may be skipped for that point.
To check for the limit of stability and for bifurcations of the response path, the tangential stiffness matrix is needed. The microplane model does not provide it directly, but it can always be computed by incrementing
the strain components (or the displacements) one by one and solving for the corresponding stress changes with the microplane model.
A greater insight into the microplane model, which may be useful for data fitting, can be achieved by separating the geometric aspect of damage (i. e., the effect of reduction of the stressresisting cross section of the material) from other inelastic phenomena. This separation has been achieved in Carol, Bazant and Prat (1991). Correlation to plasticity models and continuum damage mechanics has been elucidated in Carol and Bazant (1997).
There is another important property that is exhibited by the microplane model, and not, for example, by macroscopic plasticity models. For a nonproportional path with an abrupt change of direction such that the load increment in the ay space is directed parallel to the yield surface, the response of a plasticity model is perfectly elastic, unless this change of direction happens at a corner of the yield surface. But



Figure 14.2.3 Fitting of (a) shear compression failure envelope (in torsion) measured by Bresler and Pisler (1958), (b) biaxial failure envelope measured by Kupfer, Hilsdorf and Riisch (1969), and (c) failure envelopes in hydrostatic planes at various pressures measured by Launay and Gachon (1971). (Adapted from Bazant, Xiang el al. 1996.) 
Figure 14.2.4 Vertex effect: (a) preloading in the <тц£ц space at increasing єц and zero shear strain; (b) in the Є11Є12 space preloading corresponds to segment 01 and further tangent loading to segment 12; (c) the further tangent loading in the СГ12С12 diagram corresponds to segments 03 in classical plasticity models (fully clastic loading) and to segment 04 when vertex effect is present (after Bazant, Xiang et al. 1996). 
in reality, for all materials, this response is softer, in fact much softer, than elastic. It is as if a corner or vertex of the yield surface traveled with the state point along the path.
This effect, called the vertex effect (see Sec. 10.7 in Bazant and Cedolin 1991), is automatically described by the microplane model, but is very hard to model with the usual plastic or plasticfracturing models. It can be described only by models with many simultaneous yield surfaces, which are prohibitively difficult in the (Jij space. The microplane model is, in effect, equivalent to a set of many simultaneous yield surfaces, one for each microplane component (although these surfaces are described in the space of microplane stress components rather than in the а;у space).
This is one important advantage of the microplane approach, [t is, for example, important for obtaining the correct incremental stiffness for the case when a dc^increment (segment 12 in Fig. 14.2.4b) is superimposed on a large strain (segment 01) in the inelastic range. Segment 03 in Fig. 14.2.4c is the predicted response according to all classical macroscopic models with yield surfaces, which is elastic, and segment 04 is the prediction of microplane model, which is much softer than elastic (i. e., dcr^/dt’is < 2G where G = elastic shear modulus). Fig. 14.2.4c shows the incremental stiffness 04 calculated for the case of the present reference parameters and Єц 0.005. Indeed, the slope 04 is almost 1/5 of the slope 0.3 which would be predicted by plasticity with a simple yield surface.
The microplane model we described has been calibrated and compared to the typical test data available in the literature (Bazant, Xiang et al. 1996). They included: (1) uniaxial compression tests by van Mier (1984, 1986; Fig. 14.2.2a), for different specimen lengths and with lateral strains and volume changes measured, and by Hognestad, Hanson and McHenry (1955; Fig. 14.2.2b); (2) uniaxial direct tension tests by Petersson (1981; Fig. 14.2.2c); (3) uniaxial strain compression tests of Bazant, Bishop and Chang (1986; Fig. 14.2.2d); (4) hydrostatic compression tests by Green and Swanson (1973; Fig. 14.2.2e); (5) standard triaxial compression tests (hydrostatic loading followed by increase of one principal stress) by Balmer (1949; Fig. 14.2.2Q; (6) uniaxial cyclic compression tests of Sinha, Gerstle and Tulin (1964; Fig. 14.2.2g). (7) tests of shearcompression failure envelopes under torsion by Bresler and Pister (1958) and Goode and Helmy (1967; Fig. 14.2.3a); (8) tests of biaxial failure envelope by Kupfer, Hilsdorf and Rilsch (1969; Fig. 14.2.3b); and (9) failure envelopes from triaxial tests in octahedral plane (7rprojection) by Launay and Gachon (1971; Fig. 14.2.3c).
As seen from the figures, good fits of test data can be achieved with the microplane model. In Fig. 14.2.2a it should be noted that the uniaxial compression stressstrain diagrams are well represented for three specimens lengths, l — 5, 10, and 20 cm (it was already shown that the series coupling describes well the length effect in these tests; sec Bazant and Cedoiin 1991, Sec. 13.2). Fig. 14.2.2d serves as the basis for calibrating the volumetric stressstrain boundary, and a good fit is seen to be achieved for these enormous compressive stresses (up to 300 ksi or 2 GPa). Fig. 14.2.2f shows that the large effect of the confining pressure in standard triaxial tests can also be captured.
In Fig. 14.2.2g, note that the subsequent stress peaks in cycles reaching into the softening range are modeled quite correctly, and so are the initial unloading slopes. Significant differences, however, appear at the bottom of the cyclic loops, which is due to the fact that the unloading modulus is, in the present model, kept constant (a refinement would be possible by changing the constant unloading slope on the microplane level to a gradually decreasing slope, of course, with soine loss of simplicity). It should also be noted that the loading in these tests was quite slow and much of the curvature may have been due to relaxation caused by creep.
In Fig. 14.2.3c note that the model predicts well the shape of the failure envelopes, which is noncircular and nonhexagonal, corresponding to rounded irregular hexagons squashed from three sides. Fig. 14.2.3b shows that the ratio of uniaxial and biaxial compression strengths found in these tests can be modeled.
It must be emphasized that all the solid curves plotted in the figures are the curves that are predicted by the nticroplane model. The dotted curves in Fig. 14.2.2 arc those after correction according to the series coupling model. The dashed curves in Fig. 14.2.2c are those after correction according to the size effect law, and the dotted curves are those after a further correction according to the series coupling model.
Note that only six parameters need to be adjusted if a complete set of uniaxial, biaxial, and triaxial test data is available, and two of them can be determined separately in advance from the volumetric compression curve. If the data are limited, fewer parameters need to be adjusted. The parameters are formulated in such a manner that two of them represent scaling by affinity transformation. Normally only these two parameters need to be adjusted, which can be done by simple closedform formulas. Thus, we can conclude that the model may be efficiently used to describe concrete behavior in uniaxial, biaxial, and triaxial situations.