Mesomechanical Investigation of the Relationship between the Length of the Fracture Process Zone and Crack Extensions in Concrete

This paper deals with a mesomechanical modeling of fracture in concrete like quasi-brittle materials. The paper seeks to provide insight regarding the variation of the fracture process zone (FPZ) length during the cracking process. The correlation between the FPZ length and the crack extension is also investigated. The FPZ variation is investigated through the evolution of cohesive tangential stresses along the crack path of different notched beams under three-point bending tests. The concept of equivalent linear elastic fracture mechanics is then employed to compute the crack extensions. The mesoscale investigation shows that the relationship between the FPZ length and the crack length is non linear. Furthermore, the numerical crack extensions are used to investigate the R-curve. It is shown that the R-curve is size-dependent and notch-sensitive.


INTRODUCTION
Concrete is a typical multiscale quasi-brittle material. Research on fracture of quasi-brittle materials (concrete, rock, etc.) revealed the existence of a nonlinear zone of microcracking around the crack tip. Due to the existence of this characteristic (called the fracture process zone (FPZ)), the material behavior becomes soft and linear elastic fracture mechanics (LEFM) cannot correctly reproduce the stress field ahead the crack tip. For concrete and other quasibrittle materials, the size of the FPZ may be relatively significant compared to the specimen size. This leads to a size dependency of the strength and other fracture characteristics (fracture toughness, fracture energy, etc.). From an energetic point of view, the existence of the FPZ may be the intrinsic cause of the size effect. Hence, the features of the FPZ are important to be known for the engineering community, especially its evolution during the crack propagation.
Instead of an elastic crack tip region characterized by an infinite singular stress, a cohesive softening zone appears at the crack tip in a quasi-brittle material. Within this zone, the elastic asymptotic expansion based on LEFM is no longer sufficient. The stress distribution along the FPZ is characterized by soften-ing. Within the framework of cohesive zones models (CZM), it is assumed the existence of cohesive zone with two cohesive surfaces held together with a cohesive traction. Barenblatt [1] was the first to introduce the concept of crack-tip cohesive zone. Similar cohesive zone models have been proposed later [2]. The fictitious crack model [3] and the crack band model [4] have been introduced to study the fracture of concrete.
A good comprehension of the FPZ variation is very important to study the size effect. Both the FPZ length and the FPZ width have been a topic of many research studies both experimentally [5][6][7][8] and numerically [8][9][10]. The study of the FPZ length is essential to predict the failure of concrete members. It is particularly useful in the numerical modeling of fracture using numerical models. As far as modeling is concerned, different numerical models have been proposed in the literature. The level of description of the heterogeneities of the material is a classification criterion of the numerical models. Macroscopic models based on damage and/or plasticity theories have been extensively used to describe the fracture behavior of concrete [11,12]. At this scale, the material is considered as a continuum solid. Upon decreasing the scale, ag-gregate particles floating on a mortar matrix are observed. This so-called mesolevel (the scale of 10 -3 ) is characterized by the interactions between the aggregate particles and the cement matrix. These interactions are at the origin of local stress concentration [13]. Mesomechanical analysis of concrete behavior is now recognized as a powerful approach to predict the complex behavior of concrete (failure process [13], size effect [9,10] etc.). In a previous study [9], the authors investigated the size effect of concrete at a mesoscale level. The paper focused on the size dependency of the FPZ length. However the FPZ variation and its correlation with the crack extension during the cracking process have not been studied.
A deep investigation of the FPZ variation and its relation with the crack propagation is still needed. The correlation between the FPZ length and the crack length has been already investigated by many researchers in the literature [14,15]. However, there are no consensuses among the researchers on the evolution of the FPZ length in quasi-brittle materials. In [16,17], after the FPZ is fully open, the assumption of a constant FPZ is made. Other experimental and theoretical studies [7,18,19] show that the length of the FPZ is not always increasing during the fracture process. It increases before the FPZ is fully developed and decreases after that. So, in the scientific community, there is no widely accepted conclusion on features of variation of FPZ in concrete.
Extensive research has been conducted using cohesive zone models to study the evolution of the FPZ zone (length and width). However the relationship between the FPZ length and the crack extension needs more investigation. To study the crack extension in quasi-brittle materials, nonlinear based cohesive zone models are combined with crack propagation criteria. The stress intensity factor has been intensively used as an inherent material property [20,21]. The equivalent crack concept has been already used to study the crack extension. Irwin [22] was the first to introduce the term "equivalent crack length" to describe a fictitious increase of crack when a new distribution of stress is considered within the FPZ. The increase of the specimen compliance due to the fracture process zone development is attributed to the propagation of an equivalent elastic crack of length a eq [23] (in other words, a eq is the crack which in a specimen considered perfectly elastic, produces the same compliance as the actual specimen cracked with its damaged zone). The tip of this equivalent crack is not located at the beginning of the FPZ, but at a certain distance such that L eq = a 0 + a with a 0 is the length of the initial crack and a is the increment of the equivalent elastic crack. The crack extension is obviously related to the FPZ length. In [24], the equivalent crack length c f is supposed to be the half of the FPZ. Different concepts have been developed to study the crack propagation in quasi-brittle materials. As a consequence, different results of concrete fracture characterization are observed. As aforementioned above, there is no widely accepted conclusion about the FPZ evolution. Also, the correlation between the FPZ length and crack extension is still a subject of debate.
In the present paper, a numerical investigation of the evolution of the FPZ length in concrete using the mesoscopic approach is proposed. The objective of the present paper is twofold: (i) to provide insight into the evolution of the FPZ during the cracking process and (ii) to explore the relationship between the FPZ length and the crack extension.
The three-point bending beams tested by Rojas-Solano [25] are considered for the numerical investigations. The numerical simulations of concrete are performed at a mesoscale level using damage based model. An energetic regulation method based on the crack band approach was adopted to control the localization process [4,[26][27][28]. The exploration of the variation of the FPZ extent during the whole fracture process is based on the evolution of cohesive stress profiles obtained numerically along the crack path. The crack profiles are obtained using the procedure OUVFISS proposed by Matallah et al. [29] (implemented in the finite element software Cast3M). The crack length/FPZ length ratio is therefore investigated within the framework of equivalent LEFM. Hence, the crack extensions are used to construct the R-curves and to investigate their size dependency.

MESOSCALE INVESTIGATION OF THE GLOBAL BEHAVIOR
In the present paper, the three-point bending beams tested by Rojas-Solano [25] are considered for the numerical investigation. Mesoscale modeling of four different sizes of geometrically similar notched beams of size range 50 < D < 400 mm with notch-to-depth ratio of 0.2 (fifth notched beams) and 0.5 (half notched beams) were conducted ( Fig. 1) (more details are given in [25]).
Numerical simulation of concrete at mesoscale leads to a realistic description of the concrete behavior. The mesoscale permits an explicit representation of concrete constituents [30,31]. Concrete is considered as a biphasic material. The mortar and the ag- Fig. 1. Geometry and details of tested beams (adapted from [25]). D n = D 1 /2 n-1 , n = 1, 2, 3, 4, D 1 = 400 mm,  1 = 0.2,  2 = 0.5, b = 50 mm. gregate phases are described by their own characteristic behavior. Beyond the matrix phase, only the large aggregates are represented explicitly (the aggregates volume represents 47% of the total volume). The small aggregates and other components are assumed to be mixed up with the mortar phase establishing the matrix phase.
A softening damage law is used both for the aggregate and the mortar constituents with different characteristics. The interfacial transition zone (ITZ) is not considered. Damage creation is the result of stress concentrations occurring at the aggregate-matrix interface. The influence of the ITZ has been extensively discussed by Grondin and Matallah [30].
Numerical modeling is driven under plane stress condition with displacement control. Only the central  part of the beam where damage is expected to occur is considered with two constituents (Fig. 2). The left and right ends of the beam are considered as monophasic (macroscopic scale). Smooth transition between the different mesh zones is adopted to avoid stress concentration. A linear elastic model is considered for the macroscopic parts of the mesh to reduce the computation time. For both aggregate particles and the matrix we use an isotropic damage model proposed by Fichant [32]. The damage evolution is given by where B is a parameter which controls the slope of the softening curve defined by the exponential expression and  d0 is the strain threshold (i.e. the threshold triggering softening) related to the tensile strength ( d0 = f t /E). An energetic regulation method based on the crack band approach [4,[26][27][28] was adopted to control the mesh dependency in concrete under tension induced by the localization phenomenon. The parameter B controlling the descending branch of the softening curve is adjusted to the size of the finite element h: More details about the numerical models could be found in [9,28].

Analysis of the Global Behavior
Figure 3a and 3b show a comparison of the load-Crack Mouth Opening Displacement (CMOD) curves obtained from the experimental and the numerical simulation with the mesoscopic approach for the different beam sizes with two geometries (fifth-notched specimens and half-notched specimens) respectively. Different random draws (T1, T2 and T3) were used to generate the mesoscopic structures for each beam. An example of the crack path obtained by the postprocessing method (OUVFISS) is presented in Fig. 2.
The comparison shows that the global behavior is well reproduced numerically. The post peak behavior induced by fracture and damage are correctly reproduced for all specimen sizes. The results show nonsignificant difference of the global behavior regarding the random distribution of the aggregates because the dissipation is governed by the fracture energy. The random distribution of aggregate affects only the crack paths.

NUMERICAL INVESTIGATION OF THE FPZ LENGTH EVOLUTION
LEFM assumes the existence of an elastic crack tip region characterized by an infinite singular stress at the crack tip. Williams [33] derived the asymptotic stress field near a crack tip with the leading term exhibiting an inverse square root singularity under general planar loading conditions. Along the crack path, according to Williams series, the normal stress, which is the same as the tangential stress in polar coordinate system is given by which gives for  = 0 (mode I) the following relation: For quasi-brittle materials under mode I, a localized damage zone that contains a large number of microcracks occurs at the crack tip (the fracture process zone). Equation (2) is no longer valid to describe the behavior of the material in this zone. Instead of an elastic crack tip region characterized by an infinite singular stress, a cohesive softening zone appears at the crack tip. Within this zone, the elastic asymptotic expansion is no longer sufficient. The stress distribution along the FPZ is characterized by softening (Fig. 4). Within the framework of cohesive zones models, it is assumed the existence of cohesive zone with two cohesive surfaces held together with a cohesive traction. At a critical distance from the crack tip, the two distributions (  from LEFM and   from FPZ) are equals. To characterize the length of the FPZ, on the basis of the enrichment stress field-based criterion, the nonsingular terms of Eq. (2) have been taken into account [34,35]. The length of the FPZ is considered as the critical distance from the crack tip where the tangential stress from LEFM is equal to tensile strength f t . Based on this assumption different size effect laws have been proposed.
In the present work, the stress distribution within the FPZ region is investigated in order to follow the FPZ variation from the evolution of the tangential stress along the crack path. Numerically speaking, after the nonlinear damage computation, the crack paths are obtained using the post processing method developed by Matallah et al. [29] (Fig. 2). From the nonlinear damage finite element computation we obtain the stress values. The stress field is then expressed in the polar coordinate system to obtain the tangential stress which represents the cohesive stress. The evolution of the tangential stress   is obtained along the crack path. The FPZ starts propagation ahead the crack tip at the beginning of loading. Its length is considered as the length of the cohesive zone, i.e. the distance between the region with the cohesive stress equal to the tensile strength f t and the initial tip. When the FPZ is fully developed (i.e. the initial crack tip is stress free), the FPZ moves towards the boundary and a free stress zone of length 0 L   appears ahead the initial crack tip (see Fig. 5).

Computation of the Length of the Fully Developed FPZ
The length of the fully opened FPZ is considered as the distance between the stress-free crack tip and the zone with stress equal to the tensile strength. Figure 6 shows the evolution of tangential stress at the notch tip for each beam with respect to loading. The cohesive stress decreases progressively with further progress of cracking until it reaches zero. The FPZ length reaches its maximum value when the cohesive stress varies between zero and the tensile strength. The fully developed FPZ is formed at this time and a stress-free length occurs in front of the notch beam behind the FPZ. For each beam, the fully opened FPZ is obtained at different loading stages, i.e. the cohesive stress is reduced to zero at different loading stage. For small beams, the fully opened FPZ is formed later after the peak load. Figure 7 presents an other way of looking at the same data. The evolution of the tangential stress along the crack with respect to the height of the beam is plotted. For HN beams, the fully opened FPZ forms at 95% of the postpeak for the beam HN400, 80% of the postpeak for the beam H200, 67% of the postpeak for beam HN100 and at 22% of the postpeak for the small beams HN50. For FN beams, the fully opened FPZ forms at the peak load for the beam FN400, 90% of the postpeak for the beam F200, 70% of the postpeak for beam FN100 and at 50% of the postpeak for the small beams FN50. To derive the size effect laws, most of the theory based on LEFM assumes the development of the FPZ length at the peak load. However, as it can be seen, the fully opened FPZ corresponds to different loading stages. The plot reveals clearly that the evolution of the cohesive stress   is nonlinear along the fracture process zone. The results show notch sensitivity. For the same specimen size, the notch-to-depth ratio affects the development of the FPZ and the crack process. Figure 8 depicts the full FPZ length according to different sizes and for two initial crack lengths ( 2 = 0.2 for FN beams and ( 2 = 0.5 for HN beams). It can be seen that the full FPZ length increases with  the decrease of the relative notch length. Figure 8 affirms that the FPZ length is size-dependent and notch sensitive.

NUMERICAL ESTIMATION OF THE FPZ LENGTH VARIATION AND THE STRESS-FREE LENGTH
After the development of the fully opened FPZ, a stress-free zone occurs ahead the crack tip (Fig. 5). The crack extension depends on the FPZ length and the stress-free length. The estimation of the crack extension will be discussed in the next section.
There is a contrast among the researchers regarding the extent of the FPZ length after its fully development; if it remains constant or decreases! Numerically speaking, the best indicator to follow the FPZ size is to plot the evolution of the cohesive stress profiles along the crack path as described above. To explore the variation of the FPZ extent during the whole fracture process, the tangential stresses obtained numerically at the mesoscale level are explored at each computation step. The main objective is to give an estimation of each region indicated in Fig. 5 (the FPZ length, the stress-free length and the elastic zone).
Histograms of Figs. 9 and 10 present the formation and progress of the FPZ length L FPZ and the stress-free length 0 L    along the ligament for beams FN and HN during the entire loading process. According to this plot, the evolution of the FPZ extension can be divided into two parts: in the first part, the length of the FPZ increases gradually ones started from the notch tip until reaching its fully development whereas crack has not started yet. The second part is characterized by the crack initiation but the FPZ length is constrained by the boundary and starts to shrink. The restriction of the FPZ can be attributed to the fact that during the crack progression, the available part of the  growing. The free stress length is a part of the crack extension. The estimation of the crack length needs to investigate the portion of the FPZ that corresponds to a real crack.

ESTIMATION OF EQUIVALENT ELASTIC CRACK LENGTH
LEFM cannot be directly used for quasi-brittle materials. Nonlinear fracture behavior of concrete is due to the presence of the fracture process zone with a certain finite length. The role of the FPZ could be approximately taken into account by LEFM via the equivalent crack concept. The tip of an equivalent LEFM crack is assumed to lie ahead of the actual crack tip at a certain distance (Fig. 11). In [24] the equivalent crack length c f represents about half of the FPZ (i.e. Coef in Fig. 11 is supposed to be equal to 1/2). Other values have been proposed in the literature.
The objective of this section is to numerically determine the crack extension for notched beams within the framework of the equivalent linear elastic fracture mechanics, i.e., to find the value of the crack extension or equivalent elastic crack length which allows to find the same compliance or stiffness of the beam cracked with its FPZ length (and free stress length). After performing a nonlinear damage computation at mesoscale level, the load-displacement curve is obtained for each notched beam. Figure 12 shows a typical force-displacement curve for the beam FN200.
The stiffness of the notched beam is estimated at each point of the force-displacement curve. The second step is to simulate the same beam with different  notch-to-depth ratio >a 0 (i.e. with an equivalent crack length) until we find the crack length that makes it possible to approach the compliance or stiffness estimated in the first step. This process can be automated using an inverse analysis approach (see Fig. 13) using the MARQUARDT nonlinear algorithm [28,36]. To obtain the optimum solution, the sum of weighted squares of the errors between the measured data C(a 0 ) (the measured stiffness corresponding to each point selected in the F-D curve)) and the curve-fit value C(a eq ) is minimized. This function is called the chisquared error criterion. The inverse analysis procedure involves different steps as shown in Fig. 13. Two phases are distinguished: first, before the fully development of the FPZ, the equivalent crack length is given by L eq = a 0 + (Coef 1 L FPZ ). After the FPZ reaches its maximum length, the equivalent crack length is given by L eq = a 0 + L =0 + (Coef 2 L FPZ ) with (Coef 1 , The numerical optimization process shows that the ratio between the crack length and the FPZ length is not constant throughout the cracking process but is varying during the crack propagation process. The crack length-FPZ length relationship presents a nonlinear variation. The evolution of the equivalent crack length for HN beams (L eq = a) is plotted in Figs. 14 and 15. These plots show that FPZ length/crack length ratio is not constant. A nonlinear variation is observed. The ratio increases after the development of the fully FPZ.
A comparison of the evolution of the numerical equivalent crack length with available analytical function giving the crack length as a function of compli- ance proposed by Guinea [37] is also given in Figs. 14 and 15. Before the peak load, the analytical formula overestimates the crack length. In the postpeak behavior, a good agreement is observed between the numerical estimation and the analytical one.  (3), L eq -a 0 (equivalent LEFM) (4), L eq -a 0 (Guinea) (5) (color online). Figure 16 shows the extent of the FPZ as a function of the equivalent elastic crack length for all specimen sizes. For small beam sizes, the extent of FPZ increases linearly with the increase of the equivalent crack length exhibiting a first rising regime until reaching the maximum FPZ length corresponding to the development of the FPZ before crack starts to grow. The extent of the FPZ is followed by a decreasing regime. The FPZ drops after the crack start to increase. The reduction of the FPZ length is due to confinement of the FPZ. The crack length-ligament ratio combined with the compressive stress fields affect the development of the FPZ and restrain its free growth. This ascertainment was also observed in [19] using discrete element method. For large beams, the first regime is followed by an important longer plateau reflecting a free growth of the crack (a brittle failure). Figure 17 presents another way of looking at the same  data. The normalized (with respect to the structural ligament) FPZ length is plotted against the crack extension. Large structures present a small L FPZ /L ligament with a brittle failure, where in small beams; the FPZ could reach the entire dimension of the beams. As the ratio increases, the behavior shifts from brittle behavior (LEFM) to ductile behavior (plasticity). The effect of the relative notch depth is presented in Fig. 18. The normalized FPZ is plotted versus the normalized crack extension. For the same size, as the notch-to-depth ratio increase, the failure mode shift to brittle.

CONSTRUCTION OF R-CURVES WITH CONSIDERATION OF THE FPZ EVOLUTION
As for the LEFM, equivalent LEFM with one fracture parameter (equivalent crack length) is not able to assess the fracture process zone of concrete. The nonlinearities that take place at the crack tip while the FPZ is being formed can be predicted by using equivalent LEFM in conjunction with the variation of the fracture energy. The relationship between the fracture energy and the crack extension is called the crack-extension resistance curve (R-curve) and represents the response of the material near the crack tip to externally imposed loading. The R-curve defines the apparent increase of fracture toughness as the crack growth. When it is known, both failure load and critical crack length could be predicted. Analytical solutions have been proposed to construct the R-curve [18,23,38]. Another classical method widely used to construct the R-curve is based on the unloading compliance of load-displacement response from experimental tests [15,39,40]. For the determination of R-curve, two key parameters are needed: the fracture energy G and the crack extension a. These two parameters could be evaluated analytically or using numerical simulation. An estimation of the R-curve for the HN and FN beams described above is proposed. The crack extensions are evaluated following the process explained in the previous sections. The energy release rate is computed with an analytical expression (3).
In accordance with the linear elastic fracture mechanics, the resistance to crack growth G R (a) can be expressed from the energy release rate G(a). This energy release rate expresses the rate of energy change for each increase of small crack Ma ( [24,40]): W * is the complementary energy of the structure at a constant load. The resistance curves for HN beams and FN beams are presented in Figs. 19 and 20. The crack extensions are obtained following the numerical process exposed above. The energy restitution is estimated as a function of the dimensionless energy release rate (Eq. (3)).
The ascending part of the R-curve is well reproduced. The first upward part relating to the development of the fracture process is accompanied with a significant increase in resistance to crack growth. The second phase of the plateau value is represented only for beam 400 mm (HN400 and FN400). A typical Rcurve of a large size specimen show that for large elastic equivalent crack lengths, resistance to crack growth becomes independent of the crack length defining a plateau called the critical energy G Rc which corresponds to a critical crack length.
For the other beams (D = 50, 100, 200 mm), all the R-curve maintain rising trends without an obvious plateau. The same phenomenon has been already mentioned in [40,41]. In materials with rising R-curves without a plateau, as a crack propagates, the resistance to further crack propagation increases, and it requires a higher applied energy G in order to achieve crack extension. This tends to be the case in materials which undergo ductile fracture. As the specimen size decreases, the behavior shifts from brittle to ductile behavior. Small beams tested under three-point bending suffer from boundary conditions. Due to the confinement effect, the free development of the FPZ is limited. In [41], a beam of depth equal to 150 mm is suggested to approach the reference curve. Obviously, to obtain steady state propagation (plateau), large testing specimens geometries are required. The results also depend on the type of experimental test (bending, tension). Figure 21 shows the superposition of the R-curves (energy release rate versus the dimensionless ratio a/L ligament ) for all beam sizes. The size dependency of R-curve has already mentioned in the literature. The  results reveal the existence of a meeting point for the R-curves. Under the crack extension value corresponding to the meeting point, R-curve of large beams are above those of smaller specimens. After the meeting point, the tendency is inverted. The same phenomenon is observed for HN and FN beams. Large specimens should be tested to confirm these observations. The same phenomenon was reported by Dong et al. [20]. To evaluate the effect of the relative notch length (a 0 /D), R-curves for specimens with the same size and with two different notch-to-depth ratios are plotted in Fig. 22.
The R-curve varies with different relative crack length. For the same material and size, the energy release rate of the R-curve decreases with the increase of the initial crack length. The plots show a notch sensitivity of the R-curve. The results are also supported by experimental data [41].

CONCLUSIONS
The main goal of the present study was to investigate the evolution of the FPZ length in concrete and the relationship between the FPZ length and the crack extension. The computational modeling has been carried out at a mesoscale level using a damage-based model associated with a postprocessing method (OUVFISS). Regarding the evolution of the FPZ length, the numerical study showed that the FPZ length increases until it is fully developed and then decreases gradually after that. The FPZ evolution is affected by the boundary conditions. The second major finding was that the FPZ length/crack length ratio is not constant throughout the cracking process. The relationship between the FPZ length and the crack extension is nonlinear. Crack extensions computed by the proposed approach have been used to reproduce the R-curves. Another finding to emerge from this study is that the R-curve is size-dependent and notch-sensitive. For the cases studied, the numerical investigation revealed the existence of meeting point for R-curves of different sizes. The scope of this study was limited in terms of specimen size. Large specimens should be investigated to examine more closely and carefully the extrapolation to infinite sizes.