Calculating crack extension resistance of concrete based on a new crack propagation criterion

A crack propagation criterion was proposed for model I crack in concrete by using the initial fracture toughness K IC as an inherent material property. Based on this criterion, crack begins to propagate when the difference, between the stress intensity factors caused by the applied load KI and that by the cohesive stress KI , exceeds K ini IC . Finite element analyses was then carried out to calculate the complete load vs. crack mouth opening displacement (P-CMOD) curve, the critical crack propagation length DaC and the unstable fracture toughness K IC for notched beams under three-point bending. It was found that numerical results showed a good agreement with the experimental ones. Based on this crack propagation criterion, crack extension resistance, in terms of stress intensity factor, KR being able to consider the variation of fracture process zone (FPZ) was employed for describing crack propagation in concrete. KR is composed of K IC and K r I , which is actually equal to the driving action of crack extension. It was concluded that given the elastic modulus E, the uniaxial tensile strength ft, the fracture energy GF and K ini IC , the complete fracture process in concrete and the KR-curve of concrete can be calculated based on the numerical method. Finally, discussion was made on the effects of fracture process zone, GF and specimens geometries on KR-curve. 2012 Elsevier Ltd. Open access under CC BY license.


Introduction
Crack extension resistance curve (R-curve) has been used as the fracture criterion for decades. This approach based on R-curve was first proposed by Irwin [1] to describe the crack propagation in metals under conditions of small-scale. Since then, it has been widely applied to quasi-brittle materials like concrete. The R-curve, in terms of either strain energy release rate G or stress intensity factor (SIF) K, has already been established as a basic fracture characteristics for a given specimen geometry. Bažant and Kazemi [2,3] and Bažant and Jirásek [4] derived geometry dependent R-curve based on the size effect law observed from fracture tests. The R-curve is then obtained as the envelope of the fracture equilibrium curves for geometrically similar specimens but with various sizes. On the other hand, Ouyang et al. [5] defined the R-curve as an envelope of energy release rates with different specimen sizes but the same initial notch length. Jenq and Shah [6] proposed the two-parameter fracture model (TPFM) in which the R-curve is derived by assuming that an effective traction-free critical crack is proportional to the initial crack length. Wecharatana and Shah [7] proposed a method to calculate R-curve which is able to consider the nonlinear effect of fracture process zone length. Jenq and Shah [8] predicted the Rcurve using the constant K Ic criterion, in which K Ic was calculated based on the measured maximum load, the initial notch length using the formulas developed from linear elastic fracture mechanics (LEFM). Elices and Planas [9] and Planas et al. [10] analyzed the difference between cohesive materials and linear elastic ones using R-curve, and predicted the fracture behaviors of the cohesive 0950-0618 Ó 2012 Elsevier Ltd. http://dx.doi.org/10.1016/j.conbuildmat.2012.09.037 materials. Tvergaard and Hutchinson [11] utilized the highly refined finite element calculations and investigated the R-curve, for the crack growth in small-scale yielding under plane strain condition. Foote et al. [12] developed a theoretical model to calculate the R-curve using softening stress (r)-crack opening displacement (w) relationship in strain-softening materials.
It should be noted that the conventional R-curve can only determine the onset of crack instability. In order to describe the increase of fracture toughness during crack propagation, a new approach was proposed by Xu and Reinhardt [13] and Reinhardt and Xu [14] to evaluate the crack extension resistance in terms of considering the cohesive stress acting on the fracture process zone (FPZ). In this method, K R is obtained by combining the crack initiation fracture toughness K ini IC , which is an inherent material property, with K r I which represents the contribution of the cohesive stress along the fictitious crack zone. A semi-analytical method was applied to calculate K R -curve through using the curve of load vs. crack mouth opening displacement (P-CMOD) measured from experiment [13,14]. The result indicated that the K R -curve is an inherent material property and size-independent. However, some scholars questioned whether the K R -curve should be taken as an inherent material property. In the study of Mai [15], the SIF superposition method was applied to calculate K R -curve by assuming FPZ to be a fictitious crack, and it was concluded that K R -curve is not a material property, except for very large size specimens in which the crack-bridging zone is much smaller than the crack length and other dimensions. Furthermore, Hu and Wittmann [16] derived the extension of FPZ with the conventional compliance method in conjunction with a multi-cutting technique and proposed a new approach to calculate the K R -curve by introducing a parameter to consider the residual deformation. The results suggested that the size dependency of K R -curve is due to the difference in initial ligaments in concrete specimens. Following the studies of Xu and Reinhardt [13] and Reinhardt and Xu [14], Kumar and Barai [17,18] introduced the universal weight function to calculate the values of double-K fracture parameters. The weight function approach was also used to calculate the K R -curve [19] based on cohesive stress distribution during crack propagation. The size effect was observed from specimen with different sizes, especially during the unstable fracture stage. Furthermore, the effects of geometry and loading condition on the K R -curve in concrete were discussed by Kumar and Barai [20,21].
In the formula of K R -curves proposed by Xu and Reinhardt [13] and Reinhardt and Xu [14], K r I is a variable which is affected by FPZ length, and the relationship between fictitious crack opening displacement and stress acting on it, i.e., the cohesive law. Therefore, it is crucial to accurately determine FPZ and their associated cohesive laws for modeling the K R -curves of concrete. Moreover, it is noticed that the K R -curve is greatly affected by variation of FPZ length during crack propagation. Therefore, the difference in methods for determining FPZ length can affect the calculated results and shapes of K R -curve. Xu and Reinhardt [13] and Reinhardt and Xu [14] studied the K R -curve and FPZ length by means of fictitious crack model based on the assumption that the length of FPZ is constant when full FPZ is developed, i.e. the crack tip opening displacement (CTOD) exceeds the stress-free crack width w 0 in the softening traction-separation law. As a result, the K R -curve rises monotonically with the increase of the ratio of the effective crack length to the beam depth, i.e. a/D. However, many experimental and theoretical studies have demonstrated that the FPZ length increases before FPZ is fully developed and decreases after that. For instance, Wecharatana and Shah [7] suggested that FPZ includes both the mechanical interlock zone and the zone where microcracking and other inelastic deformation occurred. The analysis in crack propagation of cement-based composites using their method showed that the FPZ length decreases with the increase of the crack length in the loading descending branch. A more accurate analytical technique was developed by Ballarini et al. [22] to predict crack growth in cement-based composites. They also came out the same conclusion on variation of FPZ length as Wecharatana and Shah [7]. They concluded that the crack length and the displacement were so large that the interlocking effect decreased within the FPZ. Hu and Wittmann [23,24] also presented a model to predict FPZ length. Their analytical and test results both indicated that the FPZ length increases before full FPZ develops but decreases after that. Using the method proposed by Xu and Reinhardt [13] and Reinhardt and Xu [14] but considering the variation of FPZ during crack propagation, Xu et al. [25] calculated the K R -curve based on P-CMOD curve obtained from experiment. Their results also revealed that the K R -curve is size independent. It should be noted that, in their study, Xu et al. [25] assumed fracture energy G F as 100 N/m for concrete irrespective of mixtures rather than this value normally being derived from P-CMOD curve from notched beam test.
Moreover, due to sharp geometry, stress singularity may exist at the tip of crack in concrete or composite structures. For a concrete structure with cracks, crack propagation often occurs under service load. Therefore, the problem of stress singularity has been investigated in various studies [26,27], in which numerical methods are usually adopted. Several different numerical methods have been utilized to simulate crack propagation in concrete among which finite element method (FEM) is probably the most popular one. In traditional FEM, since fine meshes are needed to minimize the dependence of crack paths on initial mesh, re-meshing algorithm is necessary to accommodate crack propagation [28]. In general, re-meshing algorithm needs complex data structures to manage mesh updating as well as a well-validated crack propagation criterion. Fine crack tip meshes are also needed to represent the singularity in the stress filed, which may result in high computational cost. In order to overcome the drawbacks of re-meshing, an extended FEM [29] was developed to represent cracks by enriching the nodes in the discretization with discontinuous functions. Moreover, boundary element method (BEM) is a competitive alternative to FEM in simulating crack propagation. Combining the advantages of both FEM and BEM, a scaled boundary finite element method (SBFEM) was developed by Wolf and Song [30]. SBFEM has been successfully applied to simulate the dynamic crack propagation [31] and mixed-mode automatic crack propagation [32].
The main objective of this paper is to develop a numerical method to simulate the whole crack propagation in concrete by considering the variation of FPZ length. It was found that, using K ini IC derived from experiment as the parameter to determine crack propagation, the whole crack propagation can be simulated without the need of P-CMOD curve from experiment. On the other hand, notched concrete beams with depths D of 100, 200, and 300 mm, respectively, and initial crack ratios a 0 /D of 0.2, 0.3, 0.4, and 0.6, respectively, were tested under three-point bending and P-CMOD curves were obtained which were then compared with the numerical ones to validate the proposed numerical method. After that, the proposed method was employed to obtain K R -curve for the complete fracture process, which takes into account the cohesive stress acting on FPZ. Finally, the effect of variation of FPZ length and fracture energy on K R -curve, as well as size effect of K R -curve, was discussed.

Crack propagation criterion
In order to simulate fracture process, certain crack propagation criterion is needed which could be energy-or SIF-based. Based on the principle of energy conservation, Xie derived the energy-based cohesive crack propagation criterion [33] which states that crack propagates when the strain energy release rate exceeds the energy dissipation rate in FPZ. The criterion has been successfully utilized to simulate crack propagation in concrete [34]. On the other hand, SIF-based crack propagation criterions are also widely adopted among which the crack tip stable criterion is very popular [35]. According to this criterion, when SIF caused by deriving forces exceeds the one by cohesive forces, i.e. K I P 0, crack will propagate. It represents the competition between the crack driving forces which attempt to open the crack and the cohesive forces which attempt to close the crack. This criterion has been successfully used to simulate the crack propagation in reinforced concrete [36], mode-I and mixed-mode fracture [29,37,38] and multiple cohesive crack propagation [39] in concrete. It should be noted that the initial fracture toughness of concrete should be small enough so that it can be ignored when the criterion of K I P 0 is adopted. However, in early stage of crack propagation, crack propagation length is short so that, comparing the SIF caused by cohesive force and applied load, the initial facture toughness may not be small enough to be ignored, especially for mode-I fracture in plain concrete. Moreover, the criterion of K I P 0 can be used to determine crack propagation rather than crack initiation. Therefore, when using this criterion, crack is usually assumed to initiate when the maximum principal stress at the tip of crack exceeds uniaxial tensile strength of concrete [36][37][38][39].
Various experimental investigations have revealed that the fracture process in concrete structures included three different stages: crack initiation, stable crack propagation and unstable crack propagation or fracture. In order to differentiate the fracture process, the double-K fracture criterion was proposed which utilizes K ini IC and K un IC to describe fracture process of concrete [40,41]. Here, K ini IC is determined based on the initial cracking load, P ini , and the initial cracking length a 0 , while K un IC based on the measured maximum load, P max , and the critical effective crack length a C . However, it should be mentioned that a C is difficult to be accurately measured from experiment. In order to circumvent this difficulty in measuring a C , a linear asymptotic superposition method was introduced by Xu and Reinhardt [41] to take into account the nonlinearity of P-CMOD curve based on the following assumptions: 1. The nonlinear characteristic of the P-CMOD curve is caused by fictitious crack extension in front of a stress-free crack. 2. An effective crack consists of an equivalent elastic stress-free crack and equivalent elastic fictitious crack extension.
It is worth to mention that the linear asymptotic superposition method has been verified with both numerical and experimental results [17,18,40,41]. This method is therefore adopted in this study. Based on it, when load increases from zero to P ini , i.e., from Point O to Point A (see Fig. 1), there is no crack propagation and LEFM can be used to calculate SIF at the tip of crack. At Point A, the crack length is a 0 and SIF K P Ia at crack tip caused by the applied load P can be obtained based on P ini and a 0 . At this point, the initial fracture toughness K ini Ia is equal to K P Ia when a is equal to a 0 . When P increases from P ini to P b , i.e., from Point A to Point B (see Fig. 1), crack will propagate and a fictitious crack length 4a b is reached. Therefore, according to the linear asymptotic method, the effective crack length now becomes a b , which is equal to a 0 + 4a b . If unloading takes place at this point, load will reduce linearly down to zero, i.e., from Point B to Point O as shown in Fig. 1. After that, if reloading, load will increase linearly to P b , i.e., from Point O to Point B along Line OB as shown in Fig. 1. In other word, at Point B, the three-point bending beam can be regarded as a new virgin beam with an initial crack length a b , but it should be noted that the internal cohesive stress acting on 4a b is exactly as an equivalent externally applied load, which presents the nonlinear fracture characteristics of concrete.
Therefore, at Point B, Here, K P Ib is the SIF at the effective crack tip caused by load P b , K ini Ib is initial fracture toughness when a = a b , and K r Ib is the SIF at the effective crack tip caused by cohesive stress acting on 4a b . At Point C, the peak point of P-CMOD curve, P reaches the maximum load, P max , and the crack length a c reaches the critical crack length a C . The SIF at the effective crack tip caused by the applied load is equal to K un IC . Therefore, at Point C, At Point D, which is at the descending branch of P-CMOD curve, SIF of K P Id caused by the applied load P d is equal to K r Id , caused by cohesive stress acting on 4a d , plus the initial fracture toughness K ini Id when a = a d , i.e., Thus, based on the linear asymptotic superposition method, a complete fracture process of concrete, which is a nonlinear process, can be solved by LEFM. The P-CMOD curve of a three-point bending beam with initial crack length a 0 , under monotonous loading, can be regarded as an envelope of P-CMOD curves for a series of geometrically similar three-point bending beams with various effective crack length. According to the double K-theory [40,41], initial fracture toughness K ini IC is an inherent material property irrespective of effective crack length. Therefore, the SIF at the effective crack tip of a three-point bending notched beam can be obtained using the following simple relationship based on linear superposition principle: If load P and the corresponding crack length a are known, the crack opening displacement can be calculated using FEM, and the cohesive stress distribution can be obtained by using the cohesive stress-crack opening relationship, r-w. Then stress intensity factor K r I caused by cohesive stress can be obtained, and K P I caused by the applied load P can be directly calculated by the quarter point singular element approach. Therefore, considering the whole crack propagation process as the accumulation of many new micro-cracking initiations, for each crack length increment Da, load will reach the critical value which makes crack propagate. Each crack propagation process includes crack initiation and crack propagation. K ini IC is the index for estimating crack initiation which is the crack propagation criterion proposed in this paper. Generally, the crack propagation criterion can be described as follows: crack is in the critical state ð6Þ The above crack propagation criterion is adopted to simulate the whole fracture process in the proposed numerical method reported in this study using an iterative process which is illustrated by the program flow diagram shown in Fig. 2. The following steps are included in the iterative numerical process: 1. Calculate K ini IC based on a 0 and P ini , denoted as P(1) = P ini , from experiment. 2. Establish FEM model for three-point bending beam with crack length a i,j = a(j À 1) + Da (i = 1,2. . ., j = 2,3,. . .). Here Da is a      Table 2 Sizes, material parameters and calculated a w0 /D for L-series specimens. specified increment of crack length. i represents the load increment step during the iteration analyses in which the crack length keeps no change. j represents the crack length increment step during the iteration analyses. 3. Apply load P i,j to the notched beam and calculate cohesive stress r i,j . Calculate K P I and K r I by adjusting load P i,j until the Eq. (6) is satisfied. When K P I À K r I < K ini IC , load should be increased with a load incensement DP and cohesive stress r i,j should be recalculated based on the displacement caused by P i,j + DP. When K P I À K r I > K ini IC , load should be decreased with a load incensement DP and cohesive stress r i,j should be recalculated based on the displacement caused by P i,j À DP. Once K P I À K r I = K ini IC , the adjustment of the applied load will terminate and the adjusted load at the final step will be the actual load. 4. Repeat steps 2-3 for the next crack propagation. The iterative process terminates when the crack tip is closed to the boundary of beam or the value of the corresponding applied load becomes a negative value.
It can be seen from Fig. 2 that the applied load P(j) and the corresponding crack propagation length a(j) can be obtained in the (j)th steps. In the same step, the fracture parameters, such as K P I (j), K r I ðjÞ and CMOD(j), are also obtained. The numerical method with the iterative process was applied to analyze the B-, L-, H-, and P-series notched beams under three-point bending and their whole fracture processes were obtained. A typical example from numerical analyses is shown in Fig. 3 which illustrates the evolution of the fracture process of the L3-series specimen (see in Table 2).

Experimental verification
In order to validate of the proposed method, two series of concrete beams, with different a 0 /D (denoted as B-series) and sizes (denoted as L-series), where D is the depth of the notched beam, were tested under three-point bending and the corresponding P-CMOD curves were obtained. The set-up of three-point bending notched beam test is schematically illustrated in Fig. 4. On the other hand, P-CMOD curves were also obtained using the proposed numerical method in this study, through finite element analysis, with the crack propagation criterion mathematically described by Eqs. (5)-(7). The bilinear constitutive law for r-w, proposed by Peterson [42], was utilized to determine the cohesive stress acting on FPZ.
In this study, singular element was used to calculate SIF at the tip of crack. Because high stress gradients exist in the region around crack tip, special attention should be paid in that region. Therefore, a circle was set at the tip of crack, in which the crack tip point is the center of the circle and the crack propagation length Da is the radius of the circle. The first row of elements around the crack tip had a radius of Da/6, and the mid-side nodes were placed at the quarter points, i.e. had a radius of Da/24. In order to solve SIF at crack tip accurately, a crack propagation length Da of (1-2%) D was used in this study. For B-series specimens, the span L, depth D and width B were kept constant as 600 mm, 150 mm and 40 mm, respectively, but initial crack ratio a 0 /D varied from 0.2 to 0.6 (see Table 1). The  depth D varied from 100 mm to 300 (see Table 2). The concrete was made of aggregate with the maximal size of 10 mm and had the mechanical properties as f t = 2.3 MPa, E = 24GPa, and v = 0.2. In order to measure the initial load, four stain gauges were attached vertically in front of the precast notch on the both sides of specimen, and the distance of two strain gauges was 10 mm. When crack initiates and stars to propagate, the measured value of the strain gauge will decrease suddenly due to the fracture energy release. Therefore, the initial cracking load can be obtained according to the variation of the strain at the tip of precast notch. Initial fracture toughness can be calculated based on the initial crack length and initial load [40]. Table 1 shows the average results of different series.
The derived fracture parameters from experiments are presented in Tables 1 and 2. The complete P-CMOD curves obtained from both experiment and numerical simulation are presented in Figs. 5 and 6. In Tables 1 and 2, a w0 is the crack length corresponding to the cohesive stress when the initial crack tip decreases to zero. It can be seen that the P-CMOD curves obtained from numerical simulation agree well with those from experiment suggesting that the proposed numerical method, together with the crack propagation criterion, in this paper can be used for simulating the complete fracture process in concrete.

Crack extension curve of concrete
As aforementioned, the crack extension resistance can be expressed in terms of either strain energy release rate G or stress intensity factor K. Irwin [1] proposed an energy-based R-curve which is formulated as the following equation: where C R is the crack growth resistance, i.e. total work of fracture per unit area of crack propagation. Later on, the crack extension resistance in terms of stress intensity factor, which has been widely applied in practice, was introduced which is stated in Eq. (9).
G can be calculated from G = K 2 /E 0 according to the elastic fracture mechanics. Under plane strain and plane stress conditions, E 0 is equal to E/(1 À m 2 ) and E, respectively. Here, E is the elastic modulus and m is the Poisson's ratio.
Resistance curve, i.e., R-curve, can be presented if the stress intensity factors caused by P and Da are obtained during the crack propagation process.
About the definition of K R -curve, Xu and Reinhardt [13] and Reinhardt and Xu [14] claimed that Eq. (9) gave the stress intensity factors for different Da, but not K R -curve considering crack propagation in concrete, because material parameters relevant to concrete properties were not included in this equation. Therefore, a new K R -curve was proposed by them [13,14] in which the K R -curve was determined by combining the crack initiation toughness K ini IC which is an inherent material property of concrete, with the contribution due to the cohesive stress K r I determined by f t and a along the fictitious crack zone during the fracture process. The K R is formulated as follows: The crack propagation criterion is written as follows: Therefore, from Eqs. (10) and (11), the following can be derived: According to the crack propagation criterion, the applied load is the driving action that makes crack propagate by overcoming the resistance, i.e. the stress intensity factor caused by driving force is equal to the crack extension resistance. From Eq. (12), K R -curve can be obtained once the load and the related crack length are determined. Therefore, determining the load for different crack propagation length is a key issue, in order to obtain K R -curve, which can be achieved by using the proposed numerical method in this study combined with the proposed crack propagation criterion. For instance, the complete K R -curves of L-series specimens obtained through numerical analyses are presented in Fig. 7a-c, and the variation of FPZ with respect to specimen height is shown in Fig. 8. Here, a r is the crack length acting on the cohesive stress, i.e. FPZ length.

Effects of FPZ on K R -curve
Due to the effect of cohesive stress on fictitious crack zone, the crack extension resistance in concrete increases with the increase of the fictitious crack zone during crack propagation process. To evaluate the effects of FPZ on K R -curve, the complete fracture process in concrete are divided into four stages as shown in Fig. 7.
At the first stage, a is between a 0 and critical crack length a C . Before the applied load reaches P ini , there is no crack extension. The  FPZ length is equal to zero and K R is equal to the initial toughness. After that, a stable crack extension is first reached and the FPZ length will increase with the crack propagation. Thus, the effective crack length is equal to the sum of the initial crack length and the effective crack extension, while the crack extension resistance increases due to the increase of the FPZ length. When the applied load reaches P max , K R is equal to K un IC , and unstable crack propagation is thus initiated.
At the second stage, a is between a C and a w0 , the CTOD continues to increase up to the crack opening w 0 at zero stress, the effective crack length reaches a w0 and the cohesive stress at the initial crack tip decreases to zero. In this case, the FPZ is fully developed and its length reaches the maximum value. For L-series specimens (see Table 2), when the FPZ is fully developed, a w0 /D is equal to 0.87, 0.85 and 0.81, respectively, for L1-, L2-, and L3-series specimens. At this stage, as the crack propagates, the K R -curve ascends sharply.  Table 4 Material parameters and numerical results of P max and a w0 /D for P-series beams.  At the third stage, a w0 < a 0 < 0.93D, when the CTOD exceeds a w0 , a cohesive stress-free zone is formed in front of the initial crack tip. In this case, the effective crack consists of the initial crack, newly formed stress-free crack and FPZ. The variation of fracture process zone for L-series beams are shown in Fig. 8. Here, the sum of a 0 and newly formed stress-free crack length is denoted as a r=0 . Then, the variation of FPZ, after full FPZ is developed, becomes a key factor that affects K R value. Compared with the second stage, the rising tendency of K R -curve is not significant.
After that, there is a noteworthy phenomenon that the K Rcurves for specimens with different sizes exhibit rapid rising when a/D > 0.93. Therefore, at this fourth stage, i.e., when a/D > 0.93, the boundary effect is very significant, which is related to a/D rather than the ligament length. In the tests conducted in this research, the a/D corresponding to the initiation of boundary effect is 0.93. It can be seen from the results shown in Fig. 7 that the plateau is longer for specimens with greater size. In summary, for three-point bending beams in this test, the shape of K R -curve is affected by the FPZ variation when a < 0.93D but greatly by boundary effect when a > 0.93D. The FPZ variation during crack propagation is shown in Fig. 9.
The K R -curves for B-series specimens are presented in Fig. 10. The variation of fracture process zone for B-series beams are shown in Fig. 11. It can be seen that the plateau is not significant owing to shorter ligament length. The a w0 /D is 0.77, 0.82, 0.85 and 0.92 when a 0 /D is 0.2, 0.3, 0.4 and 0.6, respectively. Therefore, the third stage is shorter for greater a 0 /D. When a 0 /D = 0.6, a w0 /D is so closed to 0.93 that the plateau in K R -curve disappears. For this case, the K R -curve demonstrates continuous rise without a plateau.

Effect of the geometry on K R -curve
The size effects on K R -curve have been discussed by many researchers [6][7][8]. Although R-curve is defined differently by different researchers, it has been generally accepted that R-curve are both size-dependent and geometry-dependent. However, Xu and Reinhardt [13] concluded that K R -curve is independent on specimen size or the initial crack ratio, i.e., a 0 /D, based on the assumption that the FPZ length keeps constant after FPZ is fully developed. After that, using the same method but considering a varying full FPZ length, Xu et al. [25] investigated the K R -curve and got the same conclusion as Xu and Reinhardt [13]. On the other hand, Kumar and Barai [19] derived K R -curve using weight function method and found that the size effect on K R -curve is more prominent during the crack unstable propagation stage. Therefore, further study is necessary on the size effect on K R -curve due to inconsistent findings reported. The FPZ length may need to be considered as varied after it is fully developed.
In order to investigate the size effect on K R -curve, H-series three-point bending beams with initial crack ratio a 0 /D = 0.3 and the section depth varied from 200 mm to 1500 mm were analyzed by finite element analysis in which the following material parameters were adopted for concrete, f t = 3.25 MPa, G F = 100 N/m, K ini IC ¼ 0:7, and E = 30.5GPa. The a w0 /D obtained through numerical analyses are listed in Table 3.
It can be seen from Fig 12a that the K R -curves for specimens with different heights demonstrates different tendency, and curves meet approximately when a/D is about 0.57. The value of K R for smaller beams is less than those for larger beams, when a/D is less than 0.57; but it become greater for larger size beams, when a/D is greater than 0.57 which can be explained by the difference in the FPZ distribution during crack propagation in beams with different depths as shown in Fig 12b. With the same a/D, the FPZ length in smaller beams is less than that in larger ones, and K R is thus reasonably lower because of the shorter integration range at this early fracture stage when the developed FPZ length is short. With the decreasing of the beam depth, the value of a w0 /D increases. Therefore, the FPZ length is larger in beams with smaller size and K R is greater because of longer integration range at this late fracture stage when the developed FPZ is long. Moreover, the full FPZ is developed for most specimens when a/D are between 0.5 and 0.6. It may be the reason why the K R -curves meet at a point, where a/D is about 0.57.

Effects of the G F on K R -curve
The material parameters of P-series beams under three-point bending derived from experiment are listed in Table 4, among which G F of the beams are 71.4 N/m, 88.1 N/m and 95.9 N/m, respectively. Fig. 13 shows various K R , P-CMOD and FPZ curves obtained from numerical analyses for the P-series beams. With the increase of G F , the cohesive effect along FPZ is intensified. Therefore, the K R -curve clearly demonstrates the rising tendency (see Fig. 13b) and the value of a w0 /D decreases (see Table 4) correspondingly.

Conclusions
In this study, the initial fracture toughness K ini IC was adopted as the fracture propagation criterion for simulating crack propagation process in concrete. Based on this criterion, a new numerical method for calculating K R -curve, which consists of the contributions of cohesive stress along the FPZ and initiation toughness of concrete, was proposed. Subsequently, the FPZ variation was examined and its effect on K R -curve is investigated for the complete fracture process. Finally, the influence of several key fracture parameters on the K R -curve was discussed. The following conclusions can be drawn: (a) Given the elastic modulus E, the uniaxial tensile strength f t , the energy release rate G F and the initial fracture toughness K ini IC , the complete crack propagation process and the crack extension resistance K R -curve of concrete can be simulated accurately using the proposed numerical method. (b) Based on the crack propagation criterion, applied load overcomes resistance and makes crack extend during each micro-crack propagation step, i.e., the K R À Da curve is actually equated to the stress intensity factor caused by the applied load P plots related to Da. (c) The shape of K R -curve is affected by both the variation of FPZ and specimen boundary. The rising of K R -curve is significant before the FPZ reaches its full length. After that, the K R -curve will achieve a plateau with the FPZ shrinks. Ultimately, when a/D rises to around 0.93 for the beams in this research, the K R -curve exhibits a rapid rise due to boundary effect. (d) The K R -curve tends to be different for three-point bending beams with depth from 200 mm to 1500 mm, which depends on the fracture ligament length of specimen. The K R -curves meet at a point when a/D is around 0.57, and at that time the full FPZ has developed in most specimens. The value of K R for smaller size beams is less than that for larger size beams when a/D is less than 0.57, but it becomes greater for larger size beams when a/D is greater than 0.57.
Physical Sciences Research Council UK under the Grant of EP/ I031952/1 is great acknowledged.