This paper studies the e ect of ber volume fraction and ber length on fracture behaviour of reinforced concrete structures. Cementitious material is modeled using the non-ordinary state-based peridynamic theory. The non- local meshfree peridynamic method uni es mechanics of continuous and dis- continuous media and thus is well-suited in simulating discontinuities such as cracks compared to the traditional mesh-based nite element method. In addition, the non-ordinary state-based peridynamics enables us to incor- porate classical constitutive equations and damage criterion. Furthermore, in order to incorporate ber within the cementitious composite, a modi ed semi-discrete model is introduced. A two-dimensional rectangular notched model is developed to study the e ect of ber amount and ber length on mechanical response of cementitious composite and also on damage propaga- tion pattern. Sucient number of simulations are performed to statistically characterize the e ect of ber reinforcement on fracture behaviour of cemen- titious composite structures. Keywords: Peridynamics; Fiber volume fraction; Fiber length; Fracture behaviour. 1. Introduction Addition of bers improves tensile strength, ductility and load carrying capacity of cementitious materials [1, 2]. The performance of ber reinforced cementitious (FRC) structures is dependent on the ber volume fraction, orientation, and spatial distribution [3, 4]. Therefore, it is important to Corresponding author Preprint submitted to Elsevier April 15, 2019 develop practical analytical methods to investigate di erent parameters that a ect the strength of FRC structures. Studies on FRC structures are mostly limited to the nite element analysis (FEA) framework [5{7]. However, the FEA has some de ciencies which play an important role in FRC simulations. First of all, FEA fails to model the interface or bond between dissimilar materials. Also, external failure criterions are needed for crack initiation and propagation. Crack growth is generally predicted, and remeshing elements is necessary during the failure process. In recent years, nonlocal meshfree methods have been an alternative to the FEA in various continuum mechanics applications, particularly in fracture mechanics. The peridynamics, rst introduced by Silling [8], is a reformulation of the classical continuum mechanics in terms of integro-di erential equations, and thus is valid anywhere in a deformable body regardless of any discontinuity. The peridynamic method discretizes the background medium into particles and represents the material model by de ning the bond between particle points selected within a nite distance, referred to as `horizon'. The origi- nal peridynamics, known as bond-based peridynamics (BBPD), is limited to modeling materials with a Poissons ratio of 1 =4 [9]. Gerstle et al. [10] en- hanced BBPD by allowing the bond of two material points to carry not only pairwise forces, but also peridynamic moments. The resulting model, called micropolar peridynamics (MPPD), can simulate linear elastic materials with varying Poissons ratios. BBPD and MPPD have been primarily used for a fracture analysis of plain and reinforced concrete members under quasi-static [11{14] and impact loading [15, 16]. Non-ordinary state-based peridynamics (NOSBPD) is the most general- ized form of the peridynamics in which the mechanics of given point depends on the deformation state of all bonds established by surrounding points. NOSBPD can use any material model by incorporating the classical con- stitutive equations and classical damage criteria [17, 18]. Breitenfeld et al. [19] developed a NOSBPD discretization scheme for linear elastic materials and performed two-dimensional (2D) and three-dimensional (3D) stationary crack problems. Tupek et al. [20] introduced an approach that incorporates a classical damage criterion within the NOSBPD framework and applied to the simulation of a ballistic impact on aluminum panels. Yaghoobi and Chorzepa [21] proposed a meshless analysis framework based on NOSBPD to predict the response of FRC structures. NOSBPD uni es the mechanics of continuous and discontinuous media and enables the use of classical constitutive equations and classical damage criteria. In order 2 to incorporate bers within the cementitious matrix, a semi-discrete model is used in which the ber-matrix reaction is applied to the matrix material points. The ber-matrix reaction force, refer to as the ber force for short, is assumed to be equal to the ber pull-out force obtained from a micro- mechanical model. The meshless framework has been veri ed and validated through several examples. Later, Yaghoobi and Chorzepa [22] modi ed the meshless FRC modeling and studied the fracture behaviour of three point bending FRC beam with an asymmetric notch to demonstrate the e ective- ness of the proposed framework in modeling unguided and complex crack patterns. In this paper, the meshless framework introduced by Yaghoobi and Chorzepa [21] has been modi ed in order to incorporate the e ect of ber volume frac- tion, ber length and ber distribution in fracture analysis of FRC structures. A 2D model is provided with symmetric notches to investigate indistinctive damage growth path. Numerical results indicate that the proposed method- ology improves the prediction of a crack growth path in complex cementitious composite structures. Furthermore, results show that addition of bers im- proves the load carrying capacity of FRC samples by transferring loads from the stress concentration area and distributing to other areas. With the same ber volume fraction, the results indicate that greater peak load resulted from longer ber length. 2. Non-ordinary state-based peridynamics Dynamics of a given point xiin the NOSBPD formulation, depend on deformation of all bonds surrounding the point. Vector-states are introduced to describe the kinematics of the body (see Fig. 1). For example, the position vector-state X, also called the bond , from the perspective of xi, is de ned as the relative position of xiandxjin the reference con guration B0. Xhxjxii=xjxi=ij (1) The deformation vector-state, Y, represents the bond in the deformed body B. Yhxjxii=yjyi (2) where yi=xi+uiandyj=xj+uj.uiandujare deformations correspond- ing to xiandxj, respectively. Fig. 1 depicts a schematic for peridynamic 3 states. Indices indicate the correspondence of a variable to material points, Aidenotes the state Aat point xi, orAijcorresponds to the bond between xiandxj. It is also worth mentioning that in NOSBPD, the deformation vector-state, Y, takes the role of deformation gradient tensor, F, in classical local mechanics in describing kinematics of the body [21, 23]. Silling [23] introduced an approximation for evaluating the classical deformation gra- dient tensor in order to incorporate the classical constitutive equations in NOSBPD. Fi= [NiX j=1!ij(yjyi) (xjxi) dVj]:K1 i (3) where!ijdenotes an in uence function and Niis the number of points within the horizon of xi.Ka is shape tensor, expressed as follows at point xi. Ki=NiX j=1!ij(xjxi) (xjxi) dVj (4) yiyj xixjuj ui xjxixjxi xjxixjxixj xi Figure 1: Peridynamic states: Position vector-state X, deformation vector-state Y, and force vector-state T. In NOSBPD, it is assumed that points interact through long-range forces represented by force vector-state T[23]. The steady-state equilibrium equa- tions for point xiis given by summation of force vector-states over the horizon ofxi NiX j=1fTihxjxiiTjhxixjigdVj+bi= 0 (5) 4 where bis the body force and Vjis the volume of point xj. The force vector- state is obtained from Tihxjxii=!ijiK1 i(xjxi) (6) where is the rst Piola-Kirchho stress and determined through material constitutive equation. 2.1. Formulation of the force vector-state The discretization of force vector-state is provided by Breitenfeld et al. [19], and Yaghoobi and Chorzepa [21]. This study presents a modi ed for- mulation for a 2D plane stress problem. For linearly elastic materials with small deformations, the in nitesimal strain tensor = 0:5(FT+F)Ican be rewritten by de nition of peridynamic deformation gradient tensor (Eq. 3) as i= 1=2[NiX j=1!ijf(ujui) (xjxi) + (xjxi) (ujui)gdVj]:K1 i(7) The in nitesimal strain tensor in vector (or Voigt) notation takes the following form, i=K iNiX j=1!ijNijUijdVj (8) where for a 2D problem, Kis a 34 matrix and stores the arrays of inverse of shape tensor K. K i=2 4K1 11 0K1 12 0 0K12 0K1 22 K1 12K1 11K1 22K1 123 5 (9) Nis a 42 matrix which formed form arrays of position vector-state ij= [1;2] =xjxi, as show below: N=2 66410 01 20 023 775(10) 5 andUij=ujuiis the relative deformation vector. The material behavior of concrete is assumed to be linear elastic, and strain softening occurs once the concrete damages. An exponential damage model is selected to represent the strain softening and thus is used in de ning the constitutive equation [24]: = (1D)C: (11) where Cis the plain stress isotropic elastic moduli matrix, C=E (1 +)(12)2 41  0  1 0 0 0 1 =2(12)3 5 (12) Dis the damage variable ranging from zero to one, and is described as an exponential function of eq, the largest equivalent strain and determined by D=(0 : eq<0 10 eqeeq0 f0:eq0(13) where0=ft=E,ftis the tensile strength; Eis the Young's modulus; f is the parameter a ecting the slope of the softening branch. The equivalent strain is based on the modi ed von Mises de nition [25]: eq=k1 2k(12)I1+1 2ks (k1)2 (12)2I2 1+6k (1 +)2J2 (14) wherekrepresents the ratio of the compressive strength, fc, and the tensile strength,ft;denotes the Poisson's ratio, and I1andJ2are the rst and the second invariants of the strain tensor. To retrieve the term K1 i(xjxi) in Eq. 6, Qis de ned: Q ij=Q10Q2 0Q2Q1 (15) where Qij= [Q1;Q2] =K1 i(xjxi) =K1 iij. Finally, the force vector- state is obtained by Tihxjxii=!ijQ iji (16) 6 An incremental-iterative method is developed to obtain the steady-state solution for system of NOSBPD equilibrium equations. For a displace- ment loading increment, the dynamic relaxation method is adapted to it- eratively update the displacement eld. The implementation of dynamic relaxation method in solving NOSBPD equations is presented in Yaghoobi and Chorzepa [21]. 2.2. Stability Instabilities in NOSBPD analysis has been reported before in the form of unphysical propagation of damage [20, 26] and zero-energy mode [19, 27{29]. In a peridynamic body, fully damaged points are able to ow unconstrained. Therefore, consideration of the fully damaged points in computation of the peridynamic deformation gradient, F, for material points within the horizon of fully damaged regions will results in unphysically-large strains. In order to address this issue, Tupek et al. [26] modi ed the in uence function in order to surpass the e ect of fully damaged points on the neighboring points. !ij=!( ij )!D(Di;Dj) (17) where!( ij ) is the radial in uence function which weights the con- tribution from neighboring point. DiandDjare the values of the damage parameter at xiandxj, respectively, and calculated from a damage criteria, here Eq. 13. !D(Di;Dj) is used to remove the e ect of fully damage particles and is required to be zero if either xiandxjis fully damaged. Tupek et al. [26] employed speci c forms for the in uence functions as follows: !( ij ) =e(jijj )2(18) !D(Di;Dj) =1 : ifDi>DcrorDj>Dcr 0 : otherwise(19) whereDcris the critical damage and is assumed as 0.99. The zero-energy mode which is also common in other particle integration of Galerkin's weak form formulation is well de ned in [19, 27]. The addi- tion of a force-spring-like bond force into the NOSBPD force vector-state is referred in [19] as a method to control the zero-energy mode. Li et al. intro- duced a supplemented force state that eliminates the need for complicated 7 parameter adjustment in [19]. Moreover, Wu and Ren [30] introduced a stabi- lized displacement eld to surpass the zero-energy mode. This paper utilises a higher-order approximation for the nonlocal deformation gradient tensor proposed by Yaghoobi and Chorzepa [27] in order to successfully control the zero-energy mode. 3. Semi discrete ber model In semi-discrete modeling of the bers within a cementitious matrix, bers do not possess independent degrees of freedom. The ber-matrix reaction forces are applied to the background matrix mesh. Since each individual ber is implicitly represented, the semi-discrete model o ers the advantages of a fully-discrete model in determining local behavior of the bers and e ect of ber distribution. In addition, the semi-discrete model is computationally much more e ective particularly when thousands or millions of bers are present. In order to model a ber that bridges a crack, the closing forces associated with the ber-matrix reaction are transferred to the matrix. In this paper, the ber-matrix interaction or ber bridging force is governed by the ber pull-out behavior as P(), where Pis the ber pull-out force and  is the ber pull-out distance. Also, it is assumed that a ber contributes to the cementitious matrix material model when it bridges a crack, where a crack is associated with damaged matrix points. Therefore, it is necessary to contribute a ber if it crosses a damaged area. For this aim as seen in Fig. 2(a), a ber in uence area is de ned including all matrix points within a nite distance from the ber. A ber is activated only if at minimum a point within the ber in uence area is damaged (See Fig. 2(b)). In order to activate a ber, ber damage parameter, Df, is introduced. Fiber damage parameter ranges between 0 and 1 and is equal to the largest damage parameter of points around the ber. Df= 0 indicates no damage within a set of points bridged by the ber, and thus the ber force is not applied to the model. Whereas,Df>0 represents that at minimum one point around the ber is damaged. In such case, the ber is engaged in the model, and the scaled ber force, Pf, is applied to the model. Pf=DfP (20) Lumped-forces applied to the crack surface and distributed-force along the ber embedment are among viable assumptions to transfer the ber 8 bridging force to the cementitious matrix. Yaghoobi and Chorzepa [22] as- sumed that the ber force is uniformly distributed to matrix points within the ber in uence area as shown in Fig. 2(a). Details for ber force calcula- tions and the methodology to distribute the ber force to the matrix points are also given in [22]. (a) ber in undeformed body (b) ber bridging a damaged area Figure 2: Schematic for ber activation and ber force distribution. 4. Numerical results The mechanical response and fracture behavior of cementitious compos- ites is highly a ected by its matrix and ber material properties, ber type, ber amount and spatial distribution of bers as well as the distribution of initial imperfections from pre-existing cracks and pores. In addition, the complexity and inconsistency of damage growth pattern makes it dicult to study FRC composites with the traditional FEA method. On the other hand, peridynamics uni es the mechanics of continuous media, cracks, and discrete points and does not have the need for an extra criteria to nucleate cracks or to determine a crack growth path. These features, along with the nonlocal e ect of bers in FRC composites, make the proposed nonlocal approach, a powerful and convenient tool for a failure analysis of FRC composites. A 2D model (see Fig. 3) is developed to investigate the e ect of ber amount and ber length on the fracture behavior of FRC structures. In order to localize the damage, two symmetric notches have been provided in the model (2 mm in depth and 2 mmin width). The model is restrained at one side and pulled from the other side. Displacement conditions are applied to the shaded areas (5mmin width at both ends) as shown in Fig. 3. 9 Figure 3: 2D pre-notched plate under tension. The mechanical properties of cementitious matrix are highly a ected by the percentage proportions of the matrix ingredients. In this work, the me- chanical properties of cementitious matrix are taken as: E= 20000MPa , = 0:18,ft= 3MPa andfc= 30MPa . The softening parameter, f, is assumed to be 0 :01. 4.1. Mesh re nement In the rst step, a convergence study is performed to obtain an e ective discretization scheme for the plain cementitious matrix. In the convergence study, four types of peridynamic discretization scheme are used with varying spacing size, h= 2mm,h= 1mm,h= 0:66mm andh= 0:5mm. In all cases, the horizon size, , is assumed to 3 :015h. This value has been reported to ensure the nonlocality of failure at an acceptable computational expense [8]. Furthermore, in the convergence study, a modi cation on damage model is applied as follows, in order to minimize the e ect of mesh sensitivities: f=1 h(f0:5f) + 0:5f (21) Nothing that as suggested in [31], initial imperfections from pores and shrinkage are represented by reducing concrete strength of arbitrarily selected points. Thus, the tensile strength of 5% of points is reduced by 25% and another 5% are reduced by 50%. The weak points are selected arbitrary using uniform random generators. An example for weak particles map is provided in Fig. 4. For each spacing size, 60 di erent runs are performed to statistically analyze weak point distribution on load carrying capacity of plain concrete models. The load-displacement curves obtained from peridynamic simulations are provided in Fig. 5. The mean curve is also provided. As seen in Fig. 5 for all spacing sizes, the maximum load capacity (or peak force) and the softening part are slights a ected with the random selection of the weak particles. The mean curves for four di erent particle spacing size are compared in 10 Figure 4: Schematic map of weak points for h= 0:5mm, the strength of Green and Red points are reduced by 25% and 50%, respectively. Fig. 6. Furthermore, the probability function of the peak load and the variation of mean peak force by sample number are provided in Figs. 7(a) and 7(b). As seen in the gures, for smaller spacing sizes of h= 0:66mm andh= 0:50mm, results are converging. Although decreasing spacing size means increasing number of points and then computational expense, hereinafter to study the e ect of bers on damage behavior of FRC structures, the discretization with h= 0:50mmis used for accuracy. (a)h= 2:00mm (b)h= 1:00mm (c)h= 0:66mm (d)h= 0:50mm Figure 5: Load-de ection curves of plain concrete for di erent spacing sizes. To study how the ber addition a ects the performance of a FRC struc- ture, ve sets of FRC models are developed with varying number of bers and ber length as in Table 1. Noting that for each case, the correspond- ing ber volume fraction, vf, is obtained via Vf=Vt.Vf=Nflfd2 f=4 is the 11 0 0.02 0.04 0.06 0.08 0.101020304050 Displacement (mm)Force (N) h = 2.00 mm h = 1.00 mm h = 0.66 mm h = 0.50 mmFigure 6: Mean loaddisplacement curves for di erent spacing sizes. 35 37 39 41 43 4500.20.40.60.81 Peak force (N)Probability Density h = 2.00 mm h = 1.00 mm h = 0.66 mm h = 0.50 mm (a) 0 10 20 30 40 50 60353739414345 Number of runsPeak force (N) h = 2.00 mm h = 1.00 mm h = 0.66 mm h = 0.50 mm (b) Figure 7: (a) Probability function of the peak load, (b) e ect of the sample number on the mean peak load. total ber volume, where Nf,lfanddfare the number of bers, the ber length and the ber diameter, respectively. The plate volume is obtained for the plane-stress state with thickness of 1 mm. The parameters of mi- cromechanical pull-out force are adopted from [22]. Fig. 8 represents the spatial distribution of bers for FRC models in Table 1. Fibers are arbitrary positioned within the plate using a uniform distributed random numbers. As in convergence study, in order to take into account the stochastic character of weak points and ber distributions, 60 runs (for a FRC model in Table 1) are performed in which weak points and ber distributions are selected randomly. 4.2. E ect of number of bers Fig. 9 depicts the load-displacement curves for FRC1, FRC2 and FRC3. Mean load-displacement curves for FRC1, FRC2 and FRC3 are compared in Fig. 10. It is clearly illustrated that increased number of bers (or volume fraction) increases the peak load carried by the cementitious matrix. 12 FRC Model Number of bers, NfFiber length, lf(mm) Fiber volume fraction (%) FRC1 20 16 0.5 FRC2 40 16 1.0 FRC3 60 16 1.5 FRC4 53 12 1.0 FRC5 80 8 1.0 Table 1: Number of bers and ber length used in the FRC models Figure 8: Samples of ber distribution for FRC cases in Table 1 The probability function of the peak load and the residual load at dis- placement loading of 0 :1mm are provided in Figs. 11(a) and 12(a), re- spectively. It is clearly illustrated that by increasing number of bers (or volume fraction) the deviation from the mean value is increasing. Therefore, for higher ber amount, the ber distribution plays more important role to characterize the FRC models. Figs. 11(b) and 12(b) study the e ect of the simulation number on the mean peak load and mean residual load. For all three cases of FRC1, FRC2 and FRC3, the results converge when 40 or more simulations are performed. Therefore, the sample size of 60 provides statistically sucient number of simulations. 4.3. E ect of ber length In this section, for a constant value of ber volume fraction, the e ect of ber length is studied on load carrying capacity and fracture behaviour of FRC composites. Varying ber length, to achieve a constant ber volume fraction, number of bers needs to be adjusted. In addition to FRC2 ( lf= 16mm andNf= 40), two other FRC model sets of FRC4 ( lf= 12mm andNf= 53) and FRC5 ( lf= 8mm andNf= 80) are provided. The corresponded ber volume fraction for all three cases is vf= 1:0%. The load- 13 (a) FRC1,Nf= 20 (b) FRC2,Nf= 40 (c) FRC3,Nf= 60 Figure 9: Load-de ection curves of ber reinforced concrete for di erent ber numbers. 0 0.02 0.04 0.06 0.08 0.10102030405060 Displacement (mm)Force (N) No Fiber Nf = 20 Nf = 40 Nf = 60 Figure 10: Mean loaddisplacement curves for di erent ber numbers. displacement curves of 60 di erent runs are provided in Figs. (a) and (b) for FRC4 and FRC5, respectively. Furthermore, the mean load-displacement curves for FRC2, FRC4 and FRC5 are illustrated in Fig. 14. Although three sets represents the same ber volume fraction, the gure indicates that greater peak load resulted from longer ber length. The probability function of the peak load and the residual load at dis- placement loading of 0 :1mm are provided in Figs. 15(a) and 16(a), respec- tively. Figs. 15(b) and 16(b) study the e ect of the simulation number on 14 35 40 45 50 55 60 6500.050.10.150.20.250.3 Peak force (N)Probability Density Nf = 20 Nf = 40 Nf = 60(a) 0 10 20 30 40 50 604045505560 Number of runsPeak force (N) Nf = 20 Nf = 40 Nf = 60 (b) Figure 11: (a) Probability function of the peak load, (b) e ect of the sample number on the mean peak load, for di erent ber numbers. 0 5 10 15 20 25 3000.050.10.150.20.250.3 Residual force (N)Probability Density Nf = 20 Nf = 40 Nf = 60 (a) 0 10 20 30 40 50 60051015202530 Number of runsResidual force (N) Nf = 20 Nf = 40 Nf = 60 (b) Figure 12: (a) Probability function of the peak load, (b) e ect of the sample number on the mean peak load, for di erent ber numbers. (a) FRC4,lf= 12mm (b) FRC5,lf= 8mm Figure 13: Load-de ection curves of ber reinforced concrete for di erent ber numbers. the mean peak load and mean residual load. 15 0 0.02 0.04 0.06 0.08 0.10102030405060 Displacement (mm)Force (N) No Fiber lf = 16 mm lf = 12 mm lf = 8 mmFigure 14: Mean loaddisplacement curves for di erent ber numbers. 35 40 45 50 55 60 6500.050.10.150.20.250.3 Peak force (N)Probability Density lf = 16 mm lf = 12 mm lf = 8 mm (a) 0 10 20 30 40 50 604045505560 Number of runsPeak force (N) lf = 16 mm lf = 12 mm lf = 8 mm (b) Figure 15: (a) Probability function of the peak load, (b) e ect of the sample number on the mean peak load, for di erent ber numbers.. 0 5 10 15 20 25 3000.050.10.150.20.250.3 Residual force (N)Probability Density Nf = 20 Nf = 40 Nf = 60 (a) 0 10 20 30 40 50 60051015202530 Number of runsResidual force (N) Nf = 20 Nf = 40 Nf = 60 (b) Figure 16: (a) Probability function of the peak load, (b) e ect of the sample number on the mean peak load, for di erent ber numbers. 4.4. Damage patterns Fig. 17 presents a comparison between damage patterns from four cases of plain concrete, FRC1, FRC2 and FRC3. For each case, the damage pattern of four selected simulation are provided from the total of 60 simulations. 16 (a) No ber, vf= 0% (b) FRC1,vf= 0:5% (c) FRC2,vf= 1:0% (d) FRC3,vf= 1:5% Figure 17: Damage patterns for varying ber volume fraction at 0.1 mm displacement loading The gure indicates that the damage pattern for the plain concrete models includes damaged particles concentrated in the localized area between the notches. However, by increasing the number of bers (or volume fraction), the damaged area is extended to a larger area. In some samples, two discrete 17 (a) FRC2,lf= 16mm (b) FRC4,lf= 12mm (c) FRC5,lf= 8mm Figure 18: Damage patterns for varying ber length at 0.1 mm displacement loading cracks are developed or the damage is occurred away from the localized area. These results prove the fact that bers can improve the fracture behaviour of FRC composites by transferering loads from stress concentration area to other parts. The damage patterns presented in Fig. 18 are for the cases in Section 4.3 to study the e ect of ber length. The damage patterns for FRC4 ( lf= 12mm andNf= 53) resemble the damage patterns for FRC2 ( lf= 16mm andNf= 40). However, the damage patterns for FRC5 ( lf= 8mm and Nf= 80) are more concentrated to the localized area. In other word, with the same ber volume fraction, larger bers are more e ective in transferring the load from the stress concentrations, and then e ectively improve the load carrying capacity of a FRC specimen. 18 5. Conclusion The failure behavior of the ber reinforced (FRC) structures is charac- terized in the peridynamic analysis framework. A semi-discrete model is developed within the non-ordinary state-based peridynamic framework. The de nition of a nonlocal deformation gradient enables the use of the classical constitutive equations and damage models in the non-ordinary state-based peridynamic framework. Furthermore, the proposed semi-discrete model con- siders modeling the e ect of individual bers at an acceptable computational cost. The dynamic relaxation method is used to solve nonlinear peridynamic equations. The proposed approach enables to perform statistical analysis to characterize the FRC composites. The meshless peridynamic simulations of a notched FRC structure are performed to study the e ect of ber volume fraction and ber length on fracture behaviour of FRC composites. It is clearly illustrated that addition of bers improves the load carrying capacity of FRC samples by transferring loads from the stress concentration area and distributing to other areas. Furthermore, with the same ber volume frac- tion, the results indicate that greater peak load resulted from longer ber length.