Numerical Simulation on Size Effect of Fracture Toughness of Concrete Based on Mesomechanics

The fracture performance of concrete is size-dependent within a certain size range. A four-phase composite material numerical model of mesofracture considering a mortar matrix, coarse aggregates, an interfacial transition zone (ITZ) at the meso level and the initial defects of concrete was established. The initial defects were assumed to be distributed randomly in the ITZ of concrete. The numerical model of concrete mesofracture was established to simulate the fracture process of wedge splitting (WS) concrete specimens with widths of 200–2000 mm and three-point bending (3-p-b) concrete specimens with heights of 200–800 mm. The fracture process of concrete was simulated, and the peak load (Pmax) of concrete was predicted using the numerical model. Based on the simulating results, the influence of specimen size of WS and 3-p-b tests on the fracture parameters was analyzed. It was demonstrated that when the specimen size was large enough, the fracture toughness (KIC) value obtained by the linear elastic fracture mechanics formula was independent of the specimen size. Meanwhile, the improved boundary effect model (BEM) was employed to study the tensile strength (ft) and fracture toughness of concrete using the mesofracture numerical model. A discrete value of β = 1.0–1.4 was a sufficient approximation to determine the ft and KIC values of concrete.


Introduction
The mechanical behavior and fracture properties of engineering materials are significantly dimensionally dependent. As early as the 17th century, Mariotte [1] observed from strength tests of different materials that the load-bearing capacity was affected by specimen size because of the heterogeneity of materials. Griffith [2] discovered the size effects of fracture parameters of engineering materials through experiments and interpreted it as the effects of defects. A large number of results have also shown that the strength and fracture parameters of concrete obtained from experiments vary with the size of specimens [3][4][5][6][7]. Exploration of the regularity and mechanism of the size effect phenomenon has never stopped. However, by comparing a series of experimental results of specimens with various sizes, it is concluded that the size effect phenomenon exists at a certain scale range. For example, Wittmann [3] tested compact tensile specimens with section heights of 360, 720 and 1440 mm. It was concluded that the fracture energy of concrete is independent of the specimen's size when the ligament length reaches approximately 20 times the maximum aggregate diameter. It was also 2000 × 2000 mm, with thicknesses (B) of 200 mm. A displacement load and constraints were imposed on the 3-p-b and WS specimens, according to experiments, as shown in Figure 1a,b. It must be specially noted that for WS specimens, the displacement loads Fx and Fy applied in the X, Y directions were calculated by Equation (2). (c) Details of fracture process zone in two types of specimens.
(d) Details of mesostructure of concrete.  The finite element model of the concrete mesostructure was generated and fracture process was simulated by using ANSYS (14.0, PA, USA). The model parameters were determined according to the experiments given in [26], in which the coarse aggregate volume fraction was 45%, the diameters were 5-10 mm, the diameter gradation was determined by the Fuller curve and Walraven formula [27], and the locations were generated by the Monte Carlo method. Solid elements were adopted in the finite element model. Two-dimensional (2D) four-node solid elements with two degrees of freedom at each node were used. The initial defects were simulated by the solid elements with weakened properties, and the content of the initial defect was defined in Equation (1). The depth of the ITZ varied from 0.02 to 0.2 mm according to the research results [28][29][30][31][32]. The influence of ITZ depth on the simulation results of concrete strength was studied in [33] and drew the conclusion that the depth of ITZ had little influence on the failure mode and strength of concrete. Therefore, the depth of ITZ was defined to be 0.2 mm considering the simulation efficiency and research results. The mesostructure of the four-phase material model was only applied to the fracture process zone with a width of 60 mm (6 d max ), as shown in Figure 1a,b. p = N 1 /N 2 (1) where p is the initial defect content, N 1 is the number of ITZ defect units and N 2 is the total number of ITZ units.
To study the effects of size on concrete fracture performance, specimens with initial notch (a) were generated, which sizes varying from 200 to 2000 mm. Specifically, the sizes of 3-p-b specimens (W × S) were 200 × 800 mm, 300 × 1200 mm, 400 × 1600 mm, 500 × 2000 mm, 600 × 2400 mm, and 800 × 3200 mm, with thicknesses (B) of 200 mm. The sizes of WS specimens (W × L) were 200 × 200 mm, 400 × 400 mm, 600 × 600 mm, 800 × 800 mm, 1000 × 1000 mm, 1200 × 1200 mm, 1500 × 1500 mm, and 2000 × 2000 mm, with thicknesses (B) of 200 mm. A displacement load and constraints were imposed on the 3-p-b and WS specimens, according to experiments, as shown in Figure 1a,b. It must be specially noted that for Materials 2020, 13, 1370 4 of 16 WS specimens, the displacement loads F x and F y applied in the X, Y directions were calculated by Equation (2).

Constitutive Relations and Failure Criteria
At the meso-level, the coarse aggregates were regarded as a linear elastic material and a linear elastic constitutive model was adopted, as shown in Figure 2a. An elastic-brittle constitutive model was used for mortar and ITZ, as shown in Figure 2b. A nonlinear damage constitutive model (Equation (3)) was used for concrete, as shown in Figure 2c. The maximum principal stress failure criterion was used, i.e., when the maximum principal stress of the element was greater than its allowed strength, the element began to fail.
where ε is the principal strain of concrete, E is the elastic modulus, σ is the principal stress, ε f is the peak strain, ε m is the ultimate strain, σ m is the residual strength, f t is the tensile strength and n is the softening coefficient, assigned a value of one.

Constitutive Relations and Failure Criteria
At the meso-level, the coarse aggregates were regarded as a linear elastic material and a linear elastic constitutive model was adopted, as shown in Figure 2a. An elastic--brittle constitutive model was used for mortar and ITZ, as shown in Figure 2b. A nonlinear damage constitutive model (Equation (3)) was used for concrete, as shown in Figure 2c. The maximum principal stress failure criterion was used, i.e., when the maximum principal stress of the element was greater than its allowed strength, the element began to fail.
(a) Constitutive model of coarse aggregate.
(b) Constitutive model of interface and mortar.
(c) Constitutive model of concrete.
where ε is the principal strain of concrete, E is the elastic modulus, σ is the principal stress, εf is the peak strain, εm is the ultimate strain, σm is the residual strength, ft is the tensile strength and n is the softening coefficient, assigned a value of one.

Determination of Property Parameters
The mechanical properties of the primary materials are listed in Table 1. Among them, the tensile strength and elastic modulus of mortar were calculated according to empirical formulas shown in Equations (4)-(6), as described in [34,35]. The value of c/w was determined according to [26]. Considering there were difficulties in developing mechanical experiments regarding the ITZ, its properties were determined based on the research results of existing references [36].

Determination of Property Parameters
The mechanical properties of the primary materials are listed in Table 1. Among them, the tensile strength and elastic modulus of mortar were calculated according to empirical formulas shown in Equations (4)-(6), as described in [34,35]. The value of c/w was determined according to [26]. Considering there were difficulties in developing mechanical experiments regarding the ITZ, its properties were determined based on the research results of existing references [36].
where E m is the elastic modulus of mortar, f tp is the pure tensile strength of mortar, f cm is the compressive strength of mortar and c/w is the cement-water ratio of mortar.

Determination of Mesh Size
To explore the influence of element size on the simulation results of the numerical model, we analyzed the results of concrete fracture when the mesh size of the mortar and the interface were at sizes of 0.5, 1, 2, 3 and 5 mm respectively. The results were shown in Table 2 and Figure 3. It can be seen that the failure mode is not dependent on mesh size, while the peak load P max changes when mesh size varies. Combined with the results of the experiment in [26], when the mesh size was in the range from 0.5 to 2 mm, the simulating results were relatively stable and reasonable. In this paper, the mesh size of the meso structure was 1 mm, and then it gradually increased to 5 mm in the macro structure.  (6) where Em is the elastic modulus of mortar, ftp is the pure tensile strength of mortar, fcm is the compressive strength of mortar and c/w is the cement-water ratio of mortar.

Determination of Mesh Size
To explore the influence of element size on the simulation results of the numerical model, we analyzed the results of concrete fracture when the mesh size of the mortar and the interface were at sizes of 0.5, 1, 2, 3 and 5 mm respectively. The results were shown in Table 2 and Figure 3. It can be seen that the failure mode is not dependent on mesh size, while the peak load Pmax changes when mesh size varies. Combined with the results of the experiment in [26], when the mesh size was in the range from 0.5 to 2 mm, the simulating results were relatively stable and reasonable. In this paper, the mesh size of the meso structure was 1 mm, and then it gradually increased to 5 mm in the macro structure. (e) FPZ elements with 5 mm.

Cracking Process
The fracture processes of concrete 3-p-b specimens and WS specimens were simulated and compared with the experimental results of [37], as shown in Figures 4 and 5, respectively. There were no obvious differences between the fracture phenomena of the 3-p-b specimens and WS specimens of concrete for the final failure mode. When the load reached 30% of P max , the concrete specimen of 3-p-b cracked from the initial notch tip, and when the load reached 17% of P max , the concrete Materials 2020, 13, 1370 6 of 16 specimen of WS cracked from the initial notch tip-thereby, the first microcracking or strain localization occurred and the inner microcracks accumulated [38] in both types. As the load continued to increase, cracks gradually expanded and converged with the crack at the crack tip, forming a macroscopic visible main crack at the crack tip in the fracture process zone, as shown in Figures 4b and 5b [39]. Subsequently, the cracks steadily expanded as the external load gradually increases to the peak load. Then, cracks expanded unsteadily and the load began to decrease, ultimately causing the concrete fracture failure with increasing load. It could be concluded that the primary crack of WS and 3-p-b specimens propagates along the surface of coarse aggregates, so the shape, distribution, and particle size of coarse aggregates played important roles in the concrete fracture process [40]. These fracture process simulation results agreed well with the experimental results of [37,41].
Materials 2020, 13, x FOR PEER REVIEW 6 of 16 The fracture processes of concrete 3-p-b specimens and WS specimens were simulated and compared with the experimental results of [37], as shown in Figures 4 and 5, respectively. There were no obvious differences between the fracture phenomena of the 3-p-b specimens and WS specimens of concrete for the final failure mode. When the load reached 30% of Pmax, the concrete specimen of 3p-b cracked from the initial notch tip, and when the load reached 17% of Pmax, the concrete specimen of WS cracked from the initial notch tip-thereby, the first microcracking or strain localization occurred and the inner microcracks accumulated [38] in both types. As the load continued to increase, cracks gradually expanded and converged with the crack at the crack tip, forming a macroscopic visible main crack at the crack tip in the fracture process zone, as shown in Figures 4b and 5b [39]. Subsequently, the cracks steadily expanded as the external load gradually increases to the peak load. Then, cracks expanded unsteadily and the load began to decrease, ultimately causing the concrete fracture failure with increasing load. It could be concluded that the primary crack of WS and 3-p-b specimens propagates along the surface of coarse aggregates, so the shape, distribution, and particle size of coarse aggregates played important roles in the concrete fracture process [40]. These fracture process simulation results agreed well with the experimental results of [37,41].

Peak Stress
The fracture toughness and tensile strength of concrete depend on the peak load. The mesomechanical fracture model generated in Section 2 was used to simulate the concrete fracture tests of WS and 3-p-b specimens, and proportional ultimate strength fL was calculated accordingly by Equation (7) [42]. The Load-Crack Mouth Opening Displacement (P-CMOD) curves of numerical simulation results and test results of WS and 3-p-b specimens with a dimension of 200 mm in height are shown in Figure 6. Here, the values of specimen parameters, such as the ratio of the initial notch to the height of the specimen (a/W), specimen height (W), and specimen thickness (B) were calculated The fracture processes of concrete 3-p-b specimens and WS specimens were simulated and compared with the experimental results of [37], as shown in Figures 4 and 5, respectively. There were no obvious differences between the fracture phenomena of the 3-p-b specimens and WS specimens of concrete for the final failure mode. When the load reached 30% of Pmax, the concrete specimen of 3p-b cracked from the initial notch tip, and when the load reached 17% of Pmax, the concrete specimen of WS cracked from the initial notch tip-thereby, the first microcracking or strain localization occurred and the inner microcracks accumulated [38] in both types. As the load continued to increase, cracks gradually expanded and converged with the crack at the crack tip, forming a macroscopic visible main crack at the crack tip in the fracture process zone, as shown in Figures 4b and 5b [39]. Subsequently, the cracks steadily expanded as the external load gradually increases to the peak load. Then, cracks expanded unsteadily and the load began to decrease, ultimately causing the concrete fracture failure with increasing load. It could be concluded that the primary crack of WS and 3-p-b specimens propagates along the surface of coarse aggregates, so the shape, distribution, and particle size of coarse aggregates played important roles in the concrete fracture process [40]. These fracture process simulation results agreed well with the experimental results of [37,41].

Peak Stress
The fracture toughness and tensile strength of concrete depend on the peak load. The mesomechanical fracture model generated in Section 2 was used to simulate the concrete fracture tests of WS and 3-p-b specimens, and proportional ultimate strength fL was calculated accordingly by Equation (7) [42]. The Load-Crack Mouth Opening Displacement (P-CMOD) curves of numerical simulation results and test results of WS and 3-p-b specimens with a dimension of 200 mm in height are shown in Figure 6. Here, the values of specimen parameters, such as the ratio of the initial notch to the height of the specimen (a/W), specimen height (W), and specimen thickness (B) were calculated

Peak Stress
The fracture toughness and tensile strength of concrete depend on the peak load. The meso-mechanical fracture model generated in Section 2 was used to simulate the concrete fracture tests of WS and 3-p-b specimens, and proportional ultimate strength f L was calculated accordingly by Equation (7) [42]. The Load-Crack Mouth Opening Displacement (P-CMOD) curves of numerical simulation results and test results of WS and 3-p-b specimens with a dimension of 200 mm in height are shown in Figure 6. Here, the values of specimen parameters, such as the ratio of the initial notch to the height of the specimen (a/W), specimen height (W), and specimen thickness (B) were calculated according to the numerical model, which was determined using Xu's experimental data [26]. From the results presented in Figure 6 and Tables 3 and 4, it can be seen that the simulation results are more brittle. That may be due to the elastic constitutive relation of meso materials. However, the simulation results have good discreteness and agree well with the experimental results before the peak point. To analyze the influence of specimen size on proportional ultimate strength and find the independent size of concrete strength, the fracture process of concrete specimens with heights larger than 1200 mm was further simulated, as shown in Tables 3 and 4.
where h is the effective height of the specimen section, h = W − a.
results have good discreteness and agree well with the experimental results before the peak point. To analyze the influence of specimen size on proportional ultimate strength and find the independent size of concrete strength, the fracture process of concrete specimens with heights larger than 1200 mm was further simulated, as shown in Tables 3 and 4.   2   max   2 3 Bh where h is the effective height of the specimen section, h = W − a.
The failure mode and crack propagation process did not significantly change with an increase in the specimen's dimension; meanwhile, the proportional ultimate strength fL and the crack length decreased with increasing concrete specimen dimensions for specimens of both types. However, the growth of 3-p-b specimens was slower than that of WS specimens.   [43] P(kN) CMOD(mm)   The failure mode and crack propagation process did not significantly change with an increase in the specimen's dimension; meanwhile, the proportional ultimate strength f L and the crack length decreased with increasing concrete specimen dimensions for specimens of both types. However, the growth of 3-p-b specimens was slower than that of WS specimens.

Fracture Toughness
It has been concluded [4] that the fracture toughness K IC of concrete increases with increasing of specimen size, and when the specimen size is large enough, the change in K IC with size is not noticeable [44][45][46]. To verify that K IC varies with size, as well as to find the scale-independent fracture parameters of concrete, the above experimental results and numerical simulation results were substituted into Equation (8) [47,48], and the fracture toughness of concrete was calculated and is shown in Table 5. According to the results, the variation of fracture toughness K IC with specimen height W was drawn and is shown in Figure 7. It can be concluded that the fracture toughness increases with the increase in specimen size. When the height of the specimen reaches 600 mm (W/d max = 60), the trend line is basically completely within the dashed lines, the fracture toughness no longer changes significantly and is approximately 1.10 MPa·m 1/2 , independent of the size of the specimen.
where P max is the peak load, W is the beam height, B is the beam thickness and a is the initial notch length.

Discussion of the Size Effect on Fracture Toughness and Tensile Strength
It is agreed that the fracture toughness will not increase when the specimen's size is large enough [49,50]. Compared with the other classic size effect law [5,6,8,9], the boundary effect model [51][52][53][54][55] needs less input data (only peak load Pmax) to determine the size-independent tensile strength ft and fracture toughness KIC, shown as Equation (9), which is a more suitable application for the numerical model. The peak loads (Pmax) obtained from the above numerical model, the experiments of [30] and the parameters involved can be used with this equation to determine the tensile strength of concrete specimens.

( )
where a * ∞ is the characteristic crack length of material, fully determined by ft and KIC. ae is the equivalent crack length, fully determined by the specimen size, type and a, which can be described as ae = B(α)a, where B(α) is the geometric shape parameter for 3-p-b specimens [15,22,23,56,57].
Concrete is a highly heterogeneous composite material. The evolution of the concrete microfracture process zone shows that the distribution of coarse aggregates and particle size plays important roles in the crack growth. According to [15,22,23,56,57], the virtual crack length at the tip of the initial notch (Δa) at Pmax depends on the maximum size of coarse aggregates, and a discrete coefficient (β) can be introduced to describe the relationship between Δa and the maximum diameter of coarse aggregates (dmax), shown as Equation (13).
For different types of concrete specimens, the value of β can be customized. Therefore, the nominal stress can be expressed as Equation (14).

Discussion of the Size Effect on Fracture Toughness and Tensile Strength
It is agreed that the fracture toughness will not increase when the specimen's size is large enough [49,50]. Compared with the other classic size effect law [5,6,8,9], the boundary effect model [51][52][53][54][55] needs less input data (only peak load P max ) to determine the size-independent tensile strength f t and fracture toughness K IC , shown as Equation (9), which is a more suitable application for the numerical model. The peak loads (P max ) obtained from the above numerical model, the experiments of [30] and the parameters involved can be used with this equation to determine the tensile strength of concrete specimens. 1 σ 2 n (P) where a * ∞ is the characteristic crack length of material, fully determined by f t and K IC . a e is the equivalent crack length, fully determined by the specimen size, type and a, which can be described as a e = B(α)a, where B(α) is the geometric shape parameter for 3-p-b specimens [15,22,23,56,57].
Concrete is a highly heterogeneous composite material. The evolution of the concrete micro-fracture process zone shows that the distribution of coarse aggregates and particle size plays important roles in the crack growth. According to [15,22,23,56,57], the virtual crack length at the tip of the initial notch (∆a) at P max depends on the maximum size of coarse aggregates, and a discrete coefficient (β) can be introduced to describe the relationship between ∆a and the maximum diameter of coarse aggregates (d max ), shown as Equation (13).
For different types of concrete specimens, the value of β can be customized. Therefore, the nominal stress can be expressed as Equation (14). (14) where, Accordingly, the fracture toughness K IC and tensile strength f t independent of size can be regressed by measuring the peak load P max of various size specimens. The fitting results of tensile strength and fracture toughness under various conditions of ∆a are shown in Figure 8. When β ranges from 0.6 to 3.2, the tensile strength (f t ) and fracture toughness (K IC ) vary from 2.44 to 4.01 MPa and 1.06 to 1.35 MPa·m 1/2 , respectively. It is known that the ratio of concrete tensile strength to compressive strength is approximately 1/8-1/12, and the compressive strength in [26] was 29.56 MPa. It can, thus, be inferred that the value of f t varies from 2.46 to 3.70 MPa. Meanwhile, the fracture toughness of concrete (K IC ) is relatively stable at approximately 1.10 MPa·m 1/2 . It can be concluded that the value of β is 1.0-1.4, and reasonable tensile strength f t and fracture toughness K IC can be obtained for small aggregate fracture specimens. The results agree well with those of [53], which were reached based on concrete fracture experiments with various specimen types. In conclusion, the model can be used to simulate the fracture tests of small size concrete under different conditions. Based on this, the fracture toughness and tensile strength, independent of size, can be further determined by the boundary effect theory.
Materials 2020, 13, x FOR PEER REVIEW 11 of 16 For different types of concrete specimens, the value of β can be customized. Therefore, the nominal stress can be expressed as Equation (14 where, Accordingly, the fracture toughness KIC and tensile strength ft independent of size can be regressed by measuring the peak load Pmax of various size specimens. The fitting results of tensile strength and fracture toughness under various conditions of ∆a are shown in Figure 8. When β ranges from 0.6 to 3.2, the tensile strength (ft) and fracture toughness (KIC) vary from 2.44 to 4.01 MPa and 1.06 to 1.35 MPa·m 1/2 , respectively. It is known that the ratio of concrete tensile strength to compressive strength is approximately 1/8-1/12, and the compressive strength in [26] [53], which were reached based on concrete fracture experiments with various specimen types. In conclusion, the model can be used to simulate the fracture tests of small size concrete under different conditions. Based on this, the fracture toughness and tensile strength, independent of size, can be further determined by the boundary effect theory.

Conclusions
To analyze the influence of specimen size on concrete fracture parameters, a meso-mechanical fracture model of concrete considering initial defects was established, and fracture tests of concrete WS specimens and 3-p-b specimens were simulated accordingly. The tensile strength and fracture

Conclusions
To analyze the influence of specimen size on concrete fracture parameters, a meso-mechanical fracture model of concrete considering initial defects was established, and fracture tests of concrete WS specimens and 3-p-b specimens were simulated accordingly. The tensile strength and fracture toughness were determined by the boundary effect theory based on the numerical simulation results of concrete fracture tests. The main conclusions are given as follows.
(1) The numerical model of concrete mesofracture, considering initial defects, can simulate the fracture process and predict the peak load of concrete, so it is suitable for determining concrete fracture parameters and tensile strength; (2) Based on the above mesofracture numerical model, when the height of a concrete specimen reaches 600 mm (W/d max = 60), the fracture toughness K IC calculated from P max and the initial notch length according to the linear elastic fracture mechanics formula is independent of the specimen size; (3) The tensile strength (f t ) and the fracture toughness (K IC ) which are independent in specimens of concrete can be obtained by the application of the mesofracture numerical model and the BEM. This property can be well expressed by ∆a at peak load (P max ), and the relationship between ∆a and the maximum aggregate diameter (d max ) can be established by introduced a discrete coefficient (β). Discrete values of β range from 1.0 to 1.4 are a sufficient approximation to predict the f t and K IC values of concrete.