Experimental and Numerical Study of Concrete Fracture Behavior with Multiple Cracks Based on the Meso-Model

In this paper, a series of experimental and numerical studies were carried out to investigate the effect of multiple cracks on concrete fracture behavior. Seven groups of double-crack concrete three-point bending (TPB) experiments with different crack lengths and different crack distances were carried out. The experimental results showed that the bearing capacity of double-crack specimens was slightly larger than the standard specimen with one central crack. Additionally, with an increase in the second crack length or with a crack distance reduction, the concrete’s bearing capacity increased correspondingly. Based on the experiments, a numerical meso-model was developed based on applying cohesive elements. The aggregate, mortar, interface transition zone (ITZ), and potential fracture surfaces were explicitly considered in the model. In particular, cohesive elements were used to characterize the mechanical behavior of the ITZ and potential fracture surfaces. A modified constitutive concrete model was developed by considering the potential fracture surfaces’ damage relation and friction effect. The accuracy of the developed meso-model was validated through a comparison between simulation and experiments. Based on meso-models, the influence of multiple cracks on the concrete bearing capacity was investigated by analyzing the energy evolution. The analysis results showed that the bearing capacity has a linear relation with the proportion of mode II energy consumption during the fracture process, which explains why specimens with multiple cracks have a slightly larger bearing capacity than the standard specimens. In summary, this study has found that in three-point bending fracture tests primarily characterized by mode I fractures, the presence of multiple cracks near the main crack slightly enhances the load-bearing capacity of the specimens. This is attributed to a slight increase in internal energy dissipation associated with the presence of these multiple cracks.


Introduction
Concrete's fracture performance has a significant impact on the safety of bridges or other structures.When concrete cracks, it can affect the corrosion resistance of the bridge and its corresponding load-bearing capacity [1][2][3][4].Therefore, it is essential to conduct detailed research on the fracture performance of concrete.
Since the theory of fracture mechanics was first used to analyze the failure behavior of concrete, a series of research works on concrete fracture mechanics have been carried out.Researchers studied the influence of the size effect [1][2][3][4][5][6][7], the loading rate [8,9], cyclic loading [10], and corrosion [11] (and so on) on concrete fracture behavior.The fracture mechanics approach can also be used to analyze complex and composite concrete structures, such as geopolymer concrete [12], reinforced concrete [13], fiber-reinforced concrete [14], and recycled concrete [15].
In practical engineering, concrete may have multiple cracks due to the construction technology and natural factors, and these cracks often affect each other.However, there are only a few papers about the fracture behavior of concrete with multiple cracks [16].
Evaluating or analyzing the bearing capacity of concrete structures with multiple cracks is still complex.
Due to the complexity of the multiple crack problem, it is hard to analyze the fracture behavior of concrete with multiple cracks using the mathematical analysis method.Thus, the mesoscopic numerical method is more suitable for studying this problem.Up to now, many mesoscopic methods have been proposed, such as the traditional finite element method [17][18][19], the lattice model [20][21][22], the extended finite element method [23][24][25], the particle flow method [26,27], the rigid body spring method [28], and the cohesive zone model (CZM) [29][30][31][32][33][34].In particular, based on the application of cohesive elements, the CZM model is among the most advantageous numerical methods for characterizing interfacial mechanical behavior.For this reason, this kind of model can be used to simulate the fracture behavior of concrete fracture surfaces since it can accurately characterize cracks along non-prescribed trajectories.
This study conducted a series of TPB experiments to investigate the effect of multiple cracks on the mode I fracture behavior.Double-crack TPB beams with different second crack lengths and different crack distances were designed and tested.In addition, a mesomodel was developed based on the CZM and the corresponding constitutive model.Finally, based on the experimental results and the developed meso-model, the fracture process and the inner energy consumption were analyzed.

Geometry and Loading Scheme
To compare with the standard mode I fracture behavior of concrete, double-crack TPB concrete beams were designed based on previous works [35].The standard TPB beams with an initial crack length of 80 mm were chosen as the reference group, and the size of this kind of beam is length (l) × height (h) × thickness (t) = 1000 mm × 200 mm × 120 mm, as shown in Figure 1a.Based on the standard specimen, two series of experiments were designed to study the effect of second crack length (40 mm, 60 mm, 80 mm) and the crack distance (80 mm, 120 mm, 160 mm) on the fracture behavior of concrete beams; the corresponding sizes of the specimens are shown in Figure 1b,c.The design information of the experiments is listed in Table 1.In particular, TPBSTD represents the standard TPB beam, TPBSC represents the double-crack beams with different second crack lengths, and TPBCD represents the double-crack beams with different crack distances.are only a few papers about the fracture behavior of concrete with multiple cracks [16].
Evaluating or analyzing the bearing capacity of concrete structures with multiple cracks is still complex.Due to the complexity of the multiple crack problem, it is hard to analyze the fracture behavior of concrete with multiple cracks using the mathematical analysis method.Thus, the mesoscopic numerical method is more suitable for studying this problem.Up to now, many mesoscopic methods have been proposed, such as the traditional finite element method [17][18][19], the lattice model [20][21][22], the extended finite element method [23][24][25], the particle flow method [26,27], the rigid body spring method [28], and the cohesive zone model (CZM) [29][30][31][32][33][34].In particular, based on the application of cohesive elements, the CZM model is among the most advantageous numerical methods for characterizing interfacial mechanical behavior.For this reason, this kind of model can be used to simulate the fracture behavior of concrete fracture surfaces since it can accurately characterize cracks along non-prescribed trajectories.
This study conducted a series of TPB experiments to investigate the effect of multiple cracks on the mode I fracture behavior.Double-crack TPB beams with different second crack lengths and different crack distances were designed and tested.In addition, a meso-model was developed based on the CZM and the corresponding constitutive model.Finally, based on the experimental results and the developed meso-model, the fracture process and the inner energy consumption were analyzed.

Geometry and Loading Scheme
To compare with the standard mode Ⅰ fracture behavior of concrete, double-crack TPB concrete beams were designed based on previous works [35].The standard TPB beams with an initial crack length of 80 mm were chosen as the reference group, and the size of this kind of beam is length (l) × height (h) × thickness (t) = 1000 mm × 200 mm × 120 mm, as shown in Figure 1a.Based on the standard specimen, two series of experiments were designed to study the effect of second crack length (40 mm, 60 mm, 80 mm) and the crack distance (80 mm, 120 mm, 160 mm) on the fracture behavior of concrete beams; the corresponding sizes of the specimens are shown in Figure 1b,c.The design information of the experiments is listed in Table 1.In particular, TPBSTD represents the standard TPB beam, TPBSC represents the double-crack beams with different second crack lengths, and TPBCD represents the double-crack beams with different crack distances.In every experimental group, four specimens were cast to eliminate random errors, and all beams were cast at once.A steel plate with a thickness of 3 mm was used to make the initial cracks.The particle size of aggregates was distributed in the range of 5~20 mm continually.P.O 42.5 cement, sand, and water were chosen to cast the mortar matrix of the concrete.The mix proportion is cement:aggregate:sand:water = 1:1.225:2.458:0.44.Through the material property tests, the standard compressive strength, tensile strength, and elastic modulus were determined to be about 49 MPa, 3.5 MPa, and 45 GPa, respectively.
The loading scheme of the beams is shown in Figure 2. The distance between the two supports is 800 mm, and two clip-on extensometers with a 2 mm range were used to record the crack mouth opening displacement (CMOD) of the two cracks.Additionally, a force sensor with a 50 kN range was used to record the loading value.The experiments were carried out on a 500 t hydraulic testing machine through the displacement control method, and the loading rate was set to about 1.5 × 10 −3 mm/s.All measured parameters were recorded through a DH-5902 testing system (sampling frequency: 10 to 100 k Hz).In every experimental group, four specimens were cast to eliminate random errors, and all beams were cast at once.A steel plate with a thickness of 3 mm was used to make the initial cracks.The particle size of aggregates was distributed in the range of 5~20 mm continually.P.O 42.5 cement, sand, and water were chosen to cast the mortar matrix of the concrete.The mix proportion is cement:aggregate:sand:water = 1:1.225:2.458:0.44.Through the material property tests, the standard compressive strength, tensile strength, and elastic modulus were determined to be about 49 MPa, 3.5 MPa, and 45 GPa, respectively.The loading scheme of the beams is shown in Figure 2. The distance between the two supports is 800 mm, and two clip-on extensometers with a 2 mm range were used to record the crack mouth opening displacement (CMOD) of the two cracks.Additionally, a force sensor with a 50 kN range was used to record the loading value.The experiments were carried out on a 500 t hydraulic testing machine through the displacement control method, and the loading rate was set to about 1.5 × 10 −3 mm/s.All measured parameters were recorded through a DH-5902 testing system (sampling frequency: 10 to 100 k Hz).

Experimental Results
The typical fracture pattern of the TPB beams is shown in Figure 3.All the crack propagation paths of the double-crack TPB beams were almost the same as the standard TPB beams, and the second crack was not observed to propagate at the macro-scale.This result indicated that in the failure process, which is dominated by mode I fracture, the second crack hardly affected the crack propagation path of the main crack (middle crack) when the length of the second crack was less than the main crack.

Experimental Results
The typical fracture pattern of the TPB beams is shown in Figure 3.All the crack propagation paths of the double-crack TPB beams were almost the same as the standard TPB beams, and the second crack was not observed to propagate at the macro-scale.This result indicated that in the failure process, which is dominated by mode I fracture, the second crack hardly affected the crack propagation path of the main crack (middle crack) when the length of the second crack was less than the main crack.In this study, the CMOD of the main crack was called  1 , and that of the second crack was called  2 .In addition, the load value obtained from the force sensor was represented by P. The typical  Figure 5 shows the   In this study, the CMOD of the main crack was called δ 1 , and that of the second crack was called δ 2 .In addition, the load value obtained from the force sensor was represented by P. The typical P − δ 1 curves of standard TPB beams are shown in Figure 4.The curves both showed an ascending branch and a softening branch.The average peak force of these cures is 7.38 kN, and the corresponding CMOD is 0.0462 mm.In this study, the CMOD of the main crack was called  1 , and that of the second crack was called  2 .In addition, the load value obtained from the force sensor was    − curves of the corresponding dou- ble-crack TPB beams with different second crack lengths.Compared to standard beam results, the double-crack beam peak loads were slightly larger.In addition, with the increase in the second crack length, the bearing capacity of the beams increased somewhat.
As for the CMOD of the second crack,   Figure 5 shows the P − δ 1 curves and δ 2 − δ 1 curves of the corresponding double-crack TPB beams with different second crack lengths.Compared to standard beam results, the double-crack beam peak loads were slightly larger.In addition, with the increase in the second crack length, the bearing capacity of the beams increased somewhat.As for the CMOD of the second crack, δ 2 increased with the increase in the loading force at first.After the loading force reached peak value, δ 2 gradually decreased while the CMOD of the main crack δ 1 continually increased.This result indicates that the second crack only propagated in the stable fracture stage (before the load reached peak value).Additionally, with the increase in second crack length, the value of the corresponding CMOD δ 2 became larger.It should be noted that during the test, the increase in the CMOD δ 2 may be smaller than the sensitivity of the extensometer; the δ 2 − δ 1 curves may appear not as smooth as others, as shown in Figure 5b,d.To conveniently compare the influence of the second crack length and crack distance on the bearing capacity of the beams, the peak load and the corresponding CMODs were extracted and shown in Table 2.It should be noted that the group TPBSC80 is also called TPBCD80.To conveniently compare the influence of the second crack length and crack distance on the bearing capacity of the beams, the peak load and the corresponding CMODs were extracted and shown in Table 2.It should be noted that the group TPBSC80 is also called TPBCD80.

Meso-Modeling Method
The meso-modeling method proposed earlier [32,33] was adopted to analyze the fracture behavior of concrete with double parallel cracks.This method characterized the potential fracture surfaces and ITZ via zero-thickness cohesive elements.These cohesive elements were inserted into all the interfaces of solid elements, as shown in Figure 7a.There are three kinds of cohesive elements in the numerical method: (1) CE_AGG for the potential fracture surfaces of aggregate; (2) CE_MOR for the potential fracture surfaces of mortar; and (3) CE_ITZ for the ITZ of concrete.The flowchart of meso-modeling [31,32] is shown in Figure 7b; the modeling method can be divided into two steps, which are the modeling of aggregates and modeling of the cohesive zone.

Constitutive Model of Concrete Potential Fracture Surfaces (1) Single-mode damage relation
The modified constitutive model [32,33] based on the bilinear damage relation was adopted to characterize the mechanical behavior of the concrete's potential fracture surfaces.The developed constitutive model especially considered the friction effect inside the crack.Figure 8 shows the damage relations of the concrete potential fracture surfaces in the normal and tangential directions.The expression of the stresses in the normal direction can be given as follows: where t n is the normal stress; k n is the stiffness in the normal direction; δ n is the normal displacement; δ n0 is the normal damage initiation displacement; δ n f is the normal failure displacement; and D is the damage factor, which can be calculated through D =

Meso-Modeling Method
The meso-modeling method proposed earlier [32,33] was adopted to analyze the fracture behavior of concrete with double parallel cracks.This method characterized the potential fracture surfaces and ITZ via zero-thickness cohesive elements.These cohesive elements were inserted into all the interfaces of solid elements, as shown in Figure 7a.There are three kinds of cohesive elements in the numerical method: (1) CE_AGG for the potential fracture surfaces of aggregate; (2) CE_MOR for the potential fracture surfaces of mortar; and (3) CE_ITZ for the ITZ of concrete.The flowchart of meso-modeling [31,32] is shown in Figure 7b; the modeling method can be divided into two steps, which are the modeling of aggregates and modeling of the cohesive zone.

Constitutive Model of Concrete Potential Fracture Surfaces
(1) Single-mode damage relation The modified constitutive model [32,33] based on the bilinear damage relation was adopted to characterize the mechanical behavior of the concrete's potential fracture surfaces.The developed constitutive model especially considered the friction effect inside the crack.Figure 8 shows the damage relations of the concrete potential fracture surfaces in the normal and tangential directions.The expression of the stresses in the normal direction can be given as follows: where n t is the normal stress;  (2) Mixed-mode damage relation The quadric criterion [31][32][33] was adopted to define the initiation of the damage process in the mixed-mode condition, and its expression can be expressed as follows: The damage relation in the tangential direction can also be similarly given as follows: where t s is the shear stress; k s is the stiffness in the tangential direction; δ s is the tangential displacement; δ s0 is the shear damage initiation displacement; δ s f is the shear failure displacement; and D is the damage factor, which can be calculated through D = (2) Mixed-mode damage relation The quadric criterion [31][32][33] was adopted to define the initiation of the damage process in the mixed-mode condition, and its expression can be expressed as follows: where t n0 is the normal strength and t s0 is the shear strength.
The PL criterion defined the evolution of the damage, and it can be given as follows: where G n is the normal fracture energy and G s is the shear fracture energy; these two fracture energies can be calculated through the geometry relation in Figure 8. G r n is the normal energy release rate in the mixed-mode condition and G r s is the shear energy release rate in the mixed-mode condition.These four parameters can be expressed as follows: where δ r n0 is the normal relative damage initial displacement in the mixed-mode condition; δ r n f is the normal relative failure displacement in the mixed-mode state; δ r s0 is the shear relative damage initial displacement in the mixed-mode condition; and δ r s f is the shear relative failure displacement in the mixed-mode condition.
Assuming that the loading path is monotonous, by substituting ( < > are the Macaulay brackets, and this equation is suitable for the condition when δ n > 0), the geometry relation in Figure 8, t n = k n δ r n0 , and t s = k s δ r s0 into Equation (3), the relative damage initial displacements can be given as follows: Additionally, by substituting (monotonous loading path) and Equation (5) into Equation ( 4), the failure displacements can also be given as follows: Based on these calculated parameters, the total displacement δ, total initial damage displacement δ 0 , and total failure displacement δ f can be calculated as follows: Hence, the damage factor in the mixed-mode condition can be finally given as follows: (3) Friction effect When the interface is damaged, friction occurs at the interface when the interface is closed.For this reason, the friction effect should be considered precisely.In this constitutive model, the fiction stress is calculated according to the interfacial sliding condition.
(a) Interfacial not sliding In this condition, the friction stress T f can be calculated as follows: where δ slide s is the sliding displacement which has been generated and T f max is the maximum static friction stress, which can be given according to the friction law: where f is the friction coefficient.(b) Interfacial sliding In this condition, the friction stress equals the maximum static friction stress.In addition, the corresponding interfacial sliding displacement should be updated: where δ slide * s is the updated interfacial sliding displacement.
(4) Stresses in the mixed model Finally, combining the damage relation and friction effect, the stresses can be given as follows: (5) Internal energy calculation According to the stresses defined above, three corresponding kinds of energy can be extracted to analyze the internal fracture behavior, which are: (1) the normal stress work E n ; (2) the shear stress work E s ; and (3) the friction stress work E f .These energies can be calculated by: where l is the length of the zero-thickness cohesive element and t is the calculation thickness.

Numerical Analysis and Discussion of the TPB Experiments 4.1. Input Data of the Finite Element Model
The meso-model was used to analyze the concrete double-crack problem, as shown in Figure 9.In the fracture process zone, the meso-model was used to characterize the fracture behavior of concrete, while in the non-fracture zone, the macro-model was used to reflect the elastic behavior of concrete.The meso-modeling zone and the macro-modeling zone were tied together on their boundary.Through repeated trials, the length of the meso-modeling zone was set to 400 mm, and this area can completely characterize the fracture behavior of the specimen.The mesh of the numerical model is shown in Figure 10.Through the mesh independence test, the element size was set to 3 mm.One typical model (meso-modeling zone) contains about 50,000 nodes, 17,000 solid elements, and 25,000 zero-thickness cohesive elements.The mesh of the numerical model is shown in Figure 10.Through the mesh independence test, the element size was set to 3 mm.One typical model (meso-modeling zone) contains about 50,000 nodes, 17,000 solid elements, and 25,000 zero-thickness cohesive elements.
According to the experimental results, repeat trial computations, and the previous works on cohesive elements [29][30][31][32][33][34], the material parameters (meso-modeling zone) of cohesive elements were determined and are listed in Table 3, assuming the aggregate will not crack.In addition, for the meso-modeling zone, the elastic moduli of the mortar and aggregate were 40 GPa and 60 GPa, respectively, and the corresponding Poisson ratios were set to 0.22.For the macro-modeling zone, the elastic modulus and Poisson's ratio were 45 GPa and 0.22, respectively.All specimens are subjected to a concentrated displacement load at the top of the midspan, and the final displacement value was 0.5 mm.The mesh of the numerical model is shown in Figure 10.Through the mesh independence test, the element size was set to 3 mm.One typical model (meso-modeling zone) contains about 50,000 nodes, 17,000 solid elements, and 25,000 zero-thickness cohesive elements.According to the experimental results, repeat trial computations, and the previous works on cohesive elements [29][30][31][32][33][34], the material parameters (meso-modeling zone) of cohesive elements were determined and are listed in Table 3, assuming the aggregate will not crack.In addition, for the meso-modeling zone, the elastic moduli of the mortar and aggregate were 40 GPa and 60 GPa, respectively, and the corresponding Poisson ratios were set to 0.22.For the macro-modeling zone, the elastic modulus and Poisson's ratio were 45 GPa and 0.22, respectively.All specimens are subjected to a concentrated displacement load at the top of the midspan, and the final displacement value was 0.5 mm.
The numerical models were solved in the ABAQUS/EXPLICIT solver [36], with the user subroutine VUMAT in which the proposed constitutive model was implemented.The numerical models were solved in the ABAQUS/EXPLICIT solver [36], with the user subroutine VUMAT in which the proposed constitutive model was implemented.The loading time was set to 1 s to ensure a quasi-static loading condition.Three random numerical specimens were calculated in each experiment group to eliminate accidental error.Through computation, all numerical models in the same group showed a very similar result.For this reason, in each group, only one typical numerical specimen was chosen to discuss in this study.

Fracture Behavior of the Standard and Double-Crack Beams
The standard group's numerical and experimental P-δ1 curves are shown in Figure 11.The numerical results fit well with the experimental ones in terms of the curve shape and amplitude.To further investigate the fracture process and corresponding stress distribution, four typical states (pre-peak state, peak state, post-peak state, and failure state) were chosen in Figure 9, and the four states were marked by A to D.
The distribution of maximum stress during the fracture process is shown in Figure 12, and the cracks are represented by deleting the cohesive elements whose damage factor is equal to 1.The stress distribution regular pattern can be summarized as follows: (1) In the early stage of the fracture process, a stress concentration zone exists in the crack tip.(2) In the crack-propagation stage, the damaged area extends upward.Although the damaged area does not turn into a macro-crack, the stress concentration zone moves to the tip of the damaged area, which means the stress concentration zone moves further than the crack.
(3) In the final stage, the damaged area extends to the top of the specimen, and the bearing capacity of the specimen is almost lost.

Fracture Behavior of the Standard and Double-Crack Beams
The standard group's numerical and experimental P-δ1 curves are shown in Figure 11.The numerical results fit well with the experimental ones in terms of the curve shape and amplitude.To further investigate the fracture process and corresponding stress distribution, four typical states (pre-peak state, peak state, post-peak state, and failure state) were chosen in Figure 9, and the four states were marked by A to D. The distribution of maximum stress during the fracture process is shown in Figure 12, and the cracks are represented by deleting the cohesive elements whose damage factor is equal to 1.The stress distribution regular pattern can be summarized as follows: (1) In the early stage of the fracture process, a stress concentration zone exists in the crack tip.(2) In the crack-propagation stage, the damaged area extends upward.Although the damaged area does not turn into a macro-crack, the stress concentration zone moves to the tip of the damaged area, which means the stress concentration zone moves further than the crack.(3) In the final stage, the damaged area extends to the top of the specimen, and the bearing capacity of the specimen is almost lost.
It should be noted that the crack does not extend along an ideal straight line due to the random distribution of aggregate.As a result, some fracture surfaces may not separate completely.Thus, a large stress still exists in some areas without complete damage, as shown in Figure 12c,d.
The final fracture pattern of the standard beam is shown in Figure 13, and the deformation shape has been magnified 15 times.Due to the low strength of the ITZ, the   It should be noted that the crack does not extend along an ideal straight line due to the random distribution of aggregate.As a result, some fracture surfaces may not separate completely.Thus, a large stress still exists in some areas without complete damage, as shown in Figure 12c,d.
The final fracture pattern of the standard beam is shown in Figure 13, and the deformation shape has been magnified 15 times.Due to the low strength of the ITZ, the crack will preferentially pass through a nearby ITZ in the fracture process.This regular pattern leads to the randomness of the fracture propagation path, which makes the simulation results closer to the experimental ones compared to the idealized numerical macro-model.14.It can be seen that the numerical results still show an agreement with the experimental ones.As was the case with the standard group, four typical states were chosen to investigate the fracture process of the specimen.The four states are marked by A to D, as shown in Figure 14a.
The maximum stress distribution of the double-crack beam (the group TPBSC80 or TPBCD80) during the fracture process is shown in Figure 15.The corresponding regular evolution of the stress distribution can be summarized as follows: (1) In the early stage of the fracture process, there is a two-stress concentration zone at both crack tips, and the stress value of the main crack is larger than the other one.(2) In the stable crack propagation stage, both stress concentration zones extend upward, and the stress concentration zone extending from the main crack moves further than one extending from the second The typical (the group TPBSC80 or TPBCD80) numerical and experimental P − δ 1 and δ 2 − δ 1 curves of a double-crack beam are shown in Figure 14.It can be seen that the numerical results still show an agreement with the experimental ones.As was the case with the standard group, four typical states were chosen to investigate the fracture process of the specimen.The four states are marked by A to D, as shown in Figure 14a.The maximum stress distribution of the double-crack beam (the group TPBSC80 or TPBCD80) during the fracture process is shown in Figure 15.The corresponding regular evolution of the stress distribution can be summarized as follows: (1) In the early stage of the fracture process, there is a two-stress concentration zone at both crack tips, and the stress value of the main crack is larger than the other one.(2) In the stable crack propagation stage, both stress concentration zones extend upward, and the stress concentration zone extending from the main crack moves further than one extending from the second crack.
(3) In the unstable crack propagation stage (post-peak stage), the stress concentration zone of the second crack gradually vanishes, which means the second crack stops propagating and gradually closes.Meanwhile, the main crack keeps propagating.(4) In the final stage, the damaged area originating from the main crack extends to the top of the specimen, and the bearing capacity of the specimen is almost lost.
To obtain more details of the fracture process of the double-crack specimen (the group TPBSC80 or TPBCD80), the deformation of the beam at the pre-peak stage (Figure 16a) and the post-peak stage (Figure 16b) was extracted, and the damaged areas were marked in red.It can be seen that at the pre-peak stage, both CMODs open and both the initial cracks extend.In the post-peak stage, the second crack closes and the main crack continues to open.To obtain more details of the fracture process of the double-crack specimen (the group TPBSC80 or TPBCD80), the deformation of the beam at the pre-peak stage (Figure 16a) and the post-peak stage (Figure 16b) was extracted, and the damaged areas were marked in red.It can be seen that at the pre-peak stage, both CMODs open and both the

Bearing Capacity Analysis
To investigate the impact of the two cracks on the specimen's bearing capacity, first, the relation between peak force and different experiment groups was analyzed, as shown in Figure 18.The figure indicates that the peak force of the double-crack specimens is a bit larger than that of the standard specimens.Similar to the experimental results, with the increase in the second crack's length or with the decrease in the crack distance, the peak force of the beam increases correspondingly, although the increase is minimal.To reveal why multiple cracks can slightly enhance the bearing capacity of the concrete specimen, an energy analysis was carried out based on the meso-model.

Bearing Capacity Analysis
To investigate the impact of the two cracks on the specimen's bearing capacity, first, the relation between peak force and different experiment groups was analyzed, as shown in Figure 18.The figure indicates that the peak force of the double-crack specimens is a bit larger than that of the standard specimens.Similar to the experimental results, with the increase in the second crack's length or with the decrease in the crack distance, the peak force of the beam increases correspondingly, although the increase is minimal.To reveal why multiple cracks can slightly enhance the bearing capacity of the concrete specimen, an energy analysis was carried out based on the meso-model.

Bearing Capacity Analysis
To investigate the impact of the two cracks on the specimen's bearing capacity, first, the relation between peak force and different experiment groups was analyzed, as shown in Figure 18.The figure indicates that the peak force of the double-crack specimens is a bit larger than that of the standard specimens.Similar to the experimental results, with the increase in the second crack's length or with the decrease in the crack distance, the peak force of the beam increases correspondingly, although the increase is minimal.To reveal why multiple cracks can slightly enhance the bearing capacity of the concrete specimen, an energy analysis was carried out based on the meso-model.

Energy Analysis
Through the analysis, the energy consumption in the double-crack beams was very similar to the standard beam.The typical energy evolution of a standard beam is shown in Figure 19.During the fracture process, the normal stress work (mode I fracture) dominates the energy consumption, and a small amount of shear stress work (mode II fracture) also exists.In addition, friction has little effect in this study.Thus, the whole fracture process is a composite fracture process dominated by mode I fractures.To investigate the relation between the peak force and the proportion of the energy increase at the pre-peak stage, the energy increases of different groups were extracted and are listed in Table 4.

Energy Analysis
Through the analysis, the energy consumption in the double-crack beams was very similar to the standard beam.The typical energy evolution of a standard beam is shown in Figure 19.During the fracture process, the normal stress work (mode I fracture) dominates the energy consumption, and a small amount of shear stress work (mode II fracture) also exists.In addition, friction has little effect in this study.Thus, the whole fracture process is a composite fracture process dominated by mode I fractures.To investigate the relation between the peak force and the proportion of the energy increase at the pre-peak stage, the energy increases of different groups were extracted and are listed in Table 4.   Through a comparison, it can be found that as the proportion of shear stress work (mode II fractures) increases, the peak force increases correspondingly.The peak force is highly relevant to the shear stress work proportion, and a linear relationship exists between these two parameters, as shown in Figure 20.Thus, it can be inferred that due to the existence of multiple cracks, more shear stress work will be consumed during the fracture process.For concrete material, the shear strength and corresponding fracture energy are much larger than the tensile ones.For this reason, in this study, the bearing capacity of the double-crack beams was slightly higher than the standard beams.To validate the statement above, a group of numerical tests with three parallel cracks (symmetrical and asymmetrical distributions of cracks) were also carried out, and the fracture paths are shown in Figure 21.The bearing capacity of three-crack specimens (8.03 kN, 8.04 kN, respectively) is still higher than the standard one (7.84kN), regardless of whether the distribution of the crack is symmetrical or asymmetrical.

Conclusions
In this paper, a series of TPB experiments were carried out to investigate the influence of multiple cracks on concrete fracture behavior.Six groups of double-crack TPB To validate the statement above, a group of numerical tests with three parallel cracks (symmetrical and asymmetrical distributions of cracks) were also carried out, and the fracture paths are shown in Figure 21.The bearing capacity of three-crack specimens (8.03 kN, 8.04 kN, respectively) is still higher than the standard one (7.84kN), regardless of whether the distribution of the crack is symmetrical or asymmetrical.To validate the statement above, a group of numerical tests with three parallel cracks (symmetrical and asymmetrical distributions of cracks) were also carried out, and the fracture paths are shown in Figure 21.The bearing capacity of three-crack specimens (8.03 kN, 8.04 kN, respectively) is still higher than the standard one (7.84kN), regardless of whether the distribution of the crack is symmetrical or asymmetrical.

Conclusions
In this paper, a series of TPB experiments were carried out to investigate the influence of multiple cracks on concrete fracture behavior.Six groups of double-crack TPB beams were designed and tested.A numerical meso-model was developed to analyze the

Conclusions
In this paper, a series of TPB experiments were carried out to investigate the influence of multiple cracks on concrete fracture behavior.Six groups of double-crack TPB beams were designed and tested.A numerical meso-model was developed to analyze the fracture behavior of double-crack concrete specimens.In this meso-model, cohesive elements were adopted to characterize the potential fracture surfaces, and the corresponding constitutive model combining the damage relation and friction effect was developed.Through comparison, all the numerical results agreed with the experimental ones.Based on this model, the influence of multiple parallel cracks on the bearing capacity of concrete was studied by analyzing the energy consumption.The conclusions can be summarized as follows: 1.
In the mode I fracture (or composite fracture dominated by mode I fracture) condition, multiple cracks in a small zone will slightly increase the bearing capacity of the concrete.With an increase in the other crack's lengths or with a decrease in the distance between cracks, the bearing capacity increases.

2.
In terms of energy consumption, the proportion of shear stress work (mode II) is highly relevant to the bearing capacity of multiple-parallel-crack concrete.Multiple parallel cracks change the proportion of mode II fractures and finally cause an increase in the concrete bearing capacity.
Based on the conclusions above, in the mode I fracture (or composite fracture dominated by mode I fracture) condition, multiple cracks in a small zone can be equally considered as one crack from the perspective of safety design.
It should be noted that in this study, only the concrete fracture behavior with multiple parallel cracks was investigated.However, in practical engineering applications, the distribution of cracks is more random and the fracture mode is more complex.Therefore, in subsequent studies, it is necessary to further analyze the mechanical behavior of multiple cracks under a mixed-mode (mode I and mode II) fracture process.

Figure 1 .
Figure 1.Design of concrete double-crack beams: (a) standard TPB beam; (b) beams with different second crack lengths; (c) beams with different crack distances.Figure 1. Design of concrete double-crack beams: (a) standard TPB beam; (b) beams with different second crack lengths; (c) beams with different crack distances.

Figure 1 .
Figure 1.Design of concrete double-crack beams: (a) standard TPB beam; (b) beams with different second crack lengths; (c) beams with different crack distances.Figure 1. Design of concrete double-crack beams: (a) standard TPB beam; (b) beams with different second crack lengths; (c) beams with different crack distances.

Figure 2 .
Figure 2. Loading scheme of the double-crack beams: (a) designed loading scheme; (b) experimental situation.

Figure 2 .
Figure 2. Loading scheme of the double-crack beams: (a) designed loading scheme; (b) experimental situation.

Figure 3 .
Figure 3.Typical fracture behavior of a concrete specimen: (a) the typical fracture behavior of a standard specimen; (b) the typical fracture behavior of a double-crack specimen.

 − 1 P
curves of standard TPB beams are shown in Figure 4.The curves both showed an ascending branch and a softening branch.The average peak force of these cures is 7.38 kN, and the corresponding CMOD is 0.0462 mm.

Figure 3 .
Figure 3.Typical fracture behavior of a concrete specimen: (a) the typical fracture behavior of a standard specimen; (b) the typical fracture behavior of a double-crack specimen.

Figure 3 .
Figure 3.Typical fracture behavior of a concrete specimen: (a) the typical fracture behavior of a standard specimen; (b) the typical fracture behavior of a double-crack specimen.

1 P
represented by P. The typical  − curves of standard TPB beams are shown in Figure 4.The curves both showed an ascending branch and a softening branch.The average peak force of these cures is 7.38 kN, and the corresponding CMOD is 0.0462 mm.

Figure 5
Figure5shows the

Figure 7 .
Figure 7. Meso-modeling of the concrete: (a) elements of the meso-model; (b) flowchart of modeling.
n k is the stiffness in the normal direction;  n is the normal displacement;  0 n is the normal damage initiation displacement;  nf is the normal failure displacement; and D is the damage factor, which can be calculated through

Figure 8 .
Figure 8. Bilinear relation of the concrete in pure normal or shear damage conditions: (a) in the normal direction; (b) in the tangential direction.

s 22 Figure 9 .
Figure 9. Numerical model of the TPB specimens.

Figure 9 .
Figure 9. Numerical model of the TPB specimens.

Figure 9 .
Figure 9. Numerical model of the TPB specimens.

Figure 10 .
Figure 10.Mesh of the meso-modeling zone: (a) general view; (b) the zero-thickness elements of aggregates CE_AGG; (c) the zero-thickness elements of mortar CE_MOR; (d) the zero-thickness elements of ITZ CE_ITZ.

Figure 10 .
Figure 10.Mesh of the meso-modeling zone: (a) general view; (b) the zero-thickness elements of aggregates CE_AGG; (c) the zero-thickness elements of mortar CE_MOR; (d) the zero-thickness elements of ITZ CE_ITZ.

Figure 11 .
Figure 11.The  − 1 P curve of the standard group with the simulation result.

Figure 11 .Figure 12 .
Figure 11.The P − δ 1 curve of the standard group with the simulation result.

Figure 13 . 1 P
Figure 13.Fracture pattern of the standard specimen.

Figure 12 .
Figure 12.Maximum principal stress distribution in different states marked in Figure 11: (a) state A; (b) state B; (c) state C; (d) state D.

Figure 12 .
Figure 12.Maximum principal stress distribution in different states marked in Figure 11: (a) state A; (b) state B; (c) state C; (d) state D.

Figure 13 .
Figure 13.Fracture pattern of the standard specimen.

Figure 13 .
Figure 13.Fracture pattern of the standard specimen.

Figure 16 .Figure 17 . 1 P
Figure 16.Fracture pattern of a double-crack specimen (the group TPBSC80 or TPBCD80): (a) prepeak stage; (b) post-peak stage; (c) the final fracture pattern.The simulation and experiment comparisons of the other groups (the P − δ 1 curves and the δ 2 − δ 1 curves) are shown in Figure The comparison results indicate that the numerical meso-model adopted in this study can appropriately characterize the fracture behavior of concrete with multiple cracks.

Figure 18 .
Figure 18.Comparison of the peak force between different experiment groups: (a) different second crack lengths; (b) different crack distances.Figure 18.Comparison of the peak force between different experiment groups: (a) different second crack lengths; (b) different crack distances.

Figure 18 .
Figure 18.Comparison of the peak force between different experiment groups: (a) different second crack lengths; (b) different crack distances.Figure 18.Comparison of the peak force between different experiment groups: (a) different second crack lengths; (b) different crack distances.

Figure 19 .
Figure 19.Typical evolution of the internal energy (TPBSTD) in (a) the whole fracture process; (b) the pre-peak stage.

Figure 19 .
Figure 19.Typical evolution of the internal energy (TPBSTD) in (a) the whole fracture process; (b) the pre-peak stage.

Figure 20 .
Figure 20.Relation between peak force and  s E proportion.

Figure 21 .
Figure 21.The fracture pattern of three-crack specimens with (a) a symmetrical crack distribution; (b) an asymmetrical crack distribution.

Figure 20 .
Figure 20.Relation between peak force and ∆E s proportion.

Figure 20 .
Figure 20.Relation between peak force and  s E proportion.

Figure 21 .
Figure 21.The fracture pattern of three-crack specimens with (a) a symmetrical crack distribution; (b) an asymmetrical crack distribution.

Figure 21 .
Figure 21.The fracture pattern of three-crack specimens with (a) a symmetrical crack distribution; (b) an asymmetrical crack distribution.

Author Contributions:
Conceptualization, investigation, data curation, methodology, writingoriginal draft, writing-review, funding acquisition, Z.W.; writing-original draft, writing-review, formal analysis, funding acquisition, supervision, W.Z.; investigation, data curation, methodology, writing-original draft, writing-review, Y.H.All authors have read and agreed to the published version of the manuscript.Funding: This research was funded by the Natural Science Foundation of Fujian Province with the research number 2022J01931 and the Research Development Fund project of Fujian University of Technology with research number GY-Z18185.Institutional Review Board Statement: Not applicable.Informed Consent Statement: Not applicable.Displacement in the normal direction δ n0 Normal displacement at the onset of interfacial softening in the mode-I condition δ r n0 Relative normal displacement at the onset of interfacial softening in the mixed-mode condition δ n f Normal displacement at the onset of interfacial failure in the mode-I condition δ r n f Relative normal displacement at the onset of interfacial failure in the mixed-mode condition δ s Total displacement in the tangential direction δ s0 Tangential displacement at the onset of interfacial softening in the mode-II/III condition δ r s0 Relative tangential displacement at the onset of interfacial softening in the mixed-mode condition δ s f Tangential displacement at the onset of interfacial failure in the mode-II/III condition δ r s f Relative tangential displacement at the onset of interfacial failure in the mixed-mode condition δ slide s Tangential sliding displacement that has been generated during the loading process δ Total relative displacement δ 0 Total relative displacement at the onset of interfacial softening in the mixedmode condition δ f Total relative displacement at the onset of interfacial failure in the mixedmode condition TPB Three-point bending CMOD Crack mouth opening displacement CZM Cohesive zone model

Table 1 .
The geometry information of double-crack concrete TPB beams.

Table 1 .
The geometry information of double-crack concrete TPB beams.

Table 2 .
Summary of experimental results.

Table 2 .
Summary of experimental results.

Table 3 .
Material parameters applied to the concrete.

Table 4 .
Energy increments of different double cracks specimens.

Table 4 .
Energy increments of different double cracks specimens.