Inclined Fiber Pullout from a Cementitious Matrix: A Numerical Study

It is well known that fibers improve the performance of cementitious composites by acting as bridging ligaments in cracks. Such bridging behavior is often studied through fiber pullout tests. The relation between the pullout force vs. slip end displacement is characteristic of the fiber-matrix interface. However, such a relation varies significantly with the fiber inclination angle. In the current work, we establish a numerical model to simulate the entire pullout process by explicitly representing the fiber, matrix and the interface for arbitrary fiber orientations. Cohesive elements endorsed with mixed-mode fracture capacities are implemented to represent the bond-slip behavior at the interface. Contact elements with Coulomb’s friction are placed at the interface to simulate frictional contact. The bond-slip behavior is first calibrated through pull-out curves for fibers aligned with the loading direction, then validated against experimental results for steel fibers oriented at 30∘ and 60∘. Parametric studies are then performed to explore the influences of both material properties (fiber yield strength, matrix tensile strength, interfacial bond) and geometric factors (fiber diameter, embedment length and inclination angle) on the overall pullout behavior, in particular on the maximum pullout load. The proposed methodology provides the necessary pull-out curves for a fiber oriented at a given angle for multi-scale models to study fracture in fiber-reinforced cementitious materials. The novelty lies in its capacity to capture the entire pullout process for a fiber with an arbitrary inclination angle.


Introduction
Since conventional concrete tends to fail in a brittle manner under excessive loading, fibers are often added to improve its ductility and durability. Typical examples of high-performance fiber-reinforced cement-based composites are slurry infiltrated fiber concrete (SIFCON) [1][2][3], engineered cement composites (ECC) [4][5][6], steel fiber-reinforced self-compacting concrete [7,8] or high performance concrete [9], as well as ultra-high toughness cementitious composites [10]. The effectiveness of the given type of fibers, steel fibers in particular, is often assessed through a pullout test, in which the force required to pull a fiber out of the hardened concrete is measured. The transmission of this force is achieved through the interfacial bond, defined as the shear stress at the interface between the fiber and the surrounding matrix [7, [11][12][13]. In order to capture the interface failure, the stress criterion [14,15] or the energy criterion [16][17][18], as well as cohesive approaches where bond stress is determined by the relative slip between the fiber and matrix [12,16,[19][20][21] have been adopted.
Naaman et al. [12] pointed out that aligned fibers (i.e., load is applied at the fiber longitudinal direction) exhibit two types of shear bonds at the interface: elastic and frictional. When the elastic shear stress at the interface exceeds the bond strength of the interface, the bond becomes frictional in nature. This means that the bond-broken process was considered abrupt, not gradual. Based on this hypothesis, Naaman et al. [12] developed a fundamental study of the bond in fiber-reinforced cement composites, as well as the relationship between pullout curves and bond shear stress vs. slip curves. In their analysis, a complete pullout curve was analytically derived from an assumed bond-slip relationship, and the dual problem in which the bond-slip relationship was obtained from an experimental pullout curve was solved. Furthermore, they related the frictional behavior of the interface by shrink-fit and fiber-matrix misfit theory. Other examples of aligned fibers can be found in [22,23]. For instance, by focusing on the interfacial properties, Chanvillard [22] developed a micro-mechanical model to account for the different phenomena observed in a non-straight fiber. Ellis et al. [23] conducted the simulation of a single fiber taking into account the fiber morphology. Even though good agreement was obtained, only aligned fibers were considered in these two models.
When fibers are randomly distributed, however, besides bond strength and friction along the interface, additional phenomena, such as fiber bending, matrix spalling and local friction effects, are involved [24][25][26][27]. Furthermore, the contributions of these micro-mechanisms are dependent on the fiber inclination angle and the fiber material properties [13,18,[24][25][26]28,29]. The contribution of fiber bending increases with inclination angle. The fiber curvature also plays an important role in the distribution of pressure against the surrounding matrix. Due to the local curvature of the fiber and residual stress at the interface, the matrix tends to crack and spall [24,26,29]. This has a direct impact on the effective embedment length and deformation within the fiber. Moreover, snubbing friction near the fiber exit point can increase the pull-out resistance [24]. Considerable efforts have been put into modeling the mentioned phenomena involved in the pullout process of inclined fibers. Some models considered the bending mechanism of inclined fibers while ignoring the matrix spalling. For example, Mortons and Groves [30] employed an elementary beam theory to calculate the force needed to produce a plastic hinge in the fiber as it is withdrawn from the matrix. In this way, the magnitude of the additional pull-out force due to fiber inclination was accounted for. The model reproduced well the experimental observations for lower inclination angles, but failed to do so for steeper ones. Li et al. [24] assumed an exponential increase of the pullout load with fiber orientation to capture the maximum pullout load; however, detailed information of the entire pullout process was omitted. By treating the fiber as a beam on an elastic foundation with variable stiffness and the possibility of spalling, Leung and Li [18] developed a micro-mechanical model to analyze the coupled fiber bending-matrix spalling mechanism and to determine crack bridging stress in random fiber-reinforced brittle matrix composites. Afterwards, the model was extended to ductile fibers [31]. Nevertheless, friction was left out of the picture. More comprehensive models, like that of Cailleux et al. [29] and Fantilli et al. [20], involve considerable parameters, which can only be estimated through pullout tests of aligned, as well as inclined fibers. In addition, tedious numerical iterative procedures are required.
In spite of the above models, detailed simulations to consider all of the aforementioned phenomena have been seldom conducted to explicitly model the entire pullout process for fibers at an arbitrary angle and considering both fiber bending and matrix spalling. In the current work, we endeavor to do so. Cohesive models capable of representing mixed-mode fracture and Coulomb's friction are employed to simulate the interface between fiber and matrix: both the de-bonding and frictional phases. In addition, fiber bending and matrix spalling are taken into account naturally thanks to the explicit representation of the fiber, the matrix and the interface in between. The rest of the paper is organized as follows. The interface bond characteristics and matrix spalling are presented in Section 2. Model calibration and validation are given in Section 3. Numerical results and parametric studies are depicted in Sections 4 and 5. Finally, relevant conclusions are drawn in Section 6.

Bond Characterization and Matrix Spalling
In this section, two essential features that determine the pullout response of a fiber with arbitrary orientation are explained in detail; namely bond characterization and the quantification of matrix spalling. The former deals with the identification of the distinct mechanisms behind interface bond Materials 2016, 9, 800 3 of 24 deterioration and friction from pullout tests of both aligned and inclined fibers. The latter tackles the matrix failure near the fiber exit point when an inclined fiber is pulled out.

Interface Bond Characterization
The shear stress vs. slip relationship is considered to be a constitutive property of the interface. Its characterization is of primary importance for predicting both the mechanical and fracture properties of fiber-reinforced composites. Naaman et al. [11] attributed the presence and combination of distinct components to the complex nature of the bond: physical and chemical adhesion, the mechanical effect of deformed or hooked fibers, the entanglement of fibers and friction, which is greatly influenced by confinement. In the current work, we concentrate on the case of a single smooth straight fiber being pulled out from a cementitious matrix. Consequently, only the interface adhesion and friction are dealt with. By fitting with the experimental data of Leung and Shapiro [26], the contribution of internal friction to interface strength will be identified and quantified. It needs to be emphasized that the same methodology can be applied for hooked or deformed steel fibers with additional efforts to represent the detailed fiber geometry.
Either for an aligned fiber or an inclined one, there exists certain frictional resistance after debonding when being pulled out. This friction is mainly dominated by the surface roughness, and it is often assumed to be constant regardless of the fiber orientation. This component is defined as internal friction, τ f (s), which is a function of the slip displacement, s. As for inclined fibers, however, the mechanisms are quite different. Pullout load can be decomposed into a parallel force and a perpendicular one. The former pulls the fiber out, while the latter bends the fiber and changes its direction during the pullout process. In addition, due to the non-uniform compression at the interface, the brittle matrix may spall, and the fiber embedment length is reduced.
In order to take into consideration the above phenomena, we propose a generic constitutive law that has three constituents to govern the interface deterioration process as follows: where µ is Coulomb's friction coefficient, θ is the fiber inclination angle with respect to the external load direction, whereas p(θ) is the compressive stress as a function of θ. The first term, τ b (s), illustrates the interfacial bond relation due to internal physical and chemical cohesion; it is often a decaying function of the slip displacement, s. In addition, τ b (s) is cohesive in nature; thus, it can be modeled through cohesive elements [32,33] with Mode-II fracture capacity. The second term, τ f (s), is a contribution from internal friction, which resists the motion between the fiber and the matrix; it can be constant or decreasing; see Figure 1. The third term gives the shear stress due to dry friction, which plays a role only when the fiber is oriented at a non-zero angle with respect to the load direction. In the current work, the bond stress is assumed to be linear-decreasing with respect to s, whereas internal friction is assumed to be constant; namely: where τ max is the maximum shear stress resisted at the interface; it covers the effect of both the internal bond and internal friction. s b c is the critical slip displacement, when the bond is completely broken.
τ f c is the internal frictional resistance, which is a property of the interface. The constitutive law encompassed in Equations (2) and (3) highlights a gradual failure at the interface; see Figure 1. This, on the one hand, differs from the abrupt bond failure assumed by Naanman et al. [12]; on the other hand, it appears particularly important, since in most practical fiber applications, the slip displacement is less than 1 mm; see Yu et al. [34]. When working in a corrosive environment, the maximum allowed crack opening in steel fiber-reinforced composites is around

Matrix Spalling
During the pullout process of inclined fibers, due to the local curvature and stretching of the fiber segment at the free end, matrix spalling is inevitable [24,26,29]. Detailed modeling of the spalling process is not an easy task. Since the deterioration process occurs within a narrow band near the interface, excessive mesh distortion often leads to convergence problems and, thus, impedes the further modeling of the entire pullout process. Consequently, simplifications are often assumed, so that the spalling part is taken off from the numerical representation, once the matrix tensile strength is reached [27,29,35,36]. The length of the eroded matrix along the fiber direction is the so-called spalling length, denoted as L sp by Laranjeira et al. [35]. A simplified failure criterion was adopted in [35], i.e., if the spalling force imposed by fiber curvature is larger than the resistance provided by the matrix at the cracked surface, the matrix will spall. Accordingly, an estimated length of the spalled matrix can be obtained as follows: where: In the above expressions, d f stands for the fiber diameter, θ is the inclination angle, P max is the peak pullout load of an aligned fiber, whereas f t is the tensile strength of the matrix. The estimated values from Equation (4) closely match the ones measured through scanning electron microscopy (SEM) by Leung and Shapiro [26].
In the current work, Equation (4) is adopted as a first approximation to represent the matrix wedge, which is to be spalled off later on. A line of cohesive elements connects the wedge to the main part of the matrix. Trial runs are performed to determine the moment at which the matrix wedge should be deactivated. From then on, the elements within the matrix wedge cease to contribute to the overall stiffness. It bears emphasis that such a treatment is only to save computational time without detriment to the modeled phenomenon, since the first principal stresses within the matrix are checked, and the spalled length is adjusted if necessary, to ensure that the matrix would not present hyper-strength.

Finite Element Methodology: Model Calibration and Validation
As mentioned before, the objective is to explicitly model the matrix, the fiber and the interface in between. In this section, a 2D setting, as shown in Figure 2a, is chosen. Both the fiber and the matrix are represented with four-node iso-parametric elements (K-L-M-N and G-H-I-J), whereas the interface is discretized into pairs of two-node line segments (I-J and K-L). The constitutive behavior of the interface is governed by a cohesive law, as shown in Figure 2b. The initial ascending part is aimed at eliminating the numerical instability, which can be caused by near-zero slip displacements derived from machine error. Consequently, the zero-damage is set for a small slip displacement, s 0 , when the bond strength, τ c , is attained; whereas the complete damage is reached for the critical slip displacement, s c , beyond which, only frictional forces exist at the interface. Additionally implemented is a pair of contact-target elements (I-J and K-L) superposed at the same location as the interface elements, so that when the interface is under compression, no penetration should occur.
In the current work, we adopt Eq. (5) as a first approximation to represent the matrix wedge 132 which is to be spalled off later on. Then trial runs are performed to determine the moment at which the 133 matrix wedge is deactivated. From then on, the elements within the matrix wedge cease to contribute 134 to the overall stiffness. It needs to be pointed out that, such a treatment is only to save computational 135 time without the detriment to the numerical results, since the first principal stresses within the matrix 136 is checked and the spalled length is adjusted if necessary, to ensure that the matrix would not present 137 hyper-strength. 138 139 As mentioned before, we aim to explicitly model the matrix, the fibre and the interface in 140 between. In a 2D setting, as shown in Fig. 2, both the fibre and the matrix are represented with 141 four-node iso-parametric elements, whereas the interface is discretised into pairs of two-node line 142 segments. The constitutive behaviour of the interface is governed by a cohesive law, as shown in  In the current work, we adopt Eq. (5) as a first approximation to represent the matrix wedge 132 which is to be spalled off later on. Then trial runs are performed to determine the moment at which the 133 matrix wedge is deactivated. From then on, the elements within the matrix wedge cease to contribute 134 to the overall stiffness. It needs to be pointed out that, such a treatment is only to save computational 135 time without the detriment to the numerical results, since the first principal stresses within the matrix 136 is checked and the spalled length is adjusted if necessary, to ensure that the matrix would not present 137 hyper-strength. by the matrix at the cracked surface, the matrix will spall. An estimated length of the spalled matrix can be obtained according to the work of Laranjeira et al. [34] as follows:

Finite element methodology: model calibration and validation
in which, d f is the fibre diameter, q is the inclination angle, P max is the peak pullout load of the aligned 129 fibre, whereas f t is the tensile strength of the matrix. The estimated values from Eq. (5) closely match 130 the ones measured through scanning electron microscopy (SEM) by Leung and Shapiro [26].

131
In the current work, we adopt Eq. (5) as a first approximation to represent the matrix wedge 132 which is to be spalled off later on. Then trial runs are performed to determine the moment at which the 133 matrix wedge is deactivated. From then on, the elements within the matrix wedge cease to contribute 134 to the overall stiffness. It needs to be pointed out that, such a treatment is only to save computational 135 time without the detriment to the numerical results, since the first principal stresses within the matrix 136 is checked and the spalled length is adjusted if necessary, to ensure that the matrix would not present 137 hyper-strength.
in which, d f is the fibre diameter, q is the inclination angle, P max is the peak pullout load of the aligned 129 fibre, whereas f t is the tensile strength of the matrix. The estimated values from Eq. (5) closely match 130 the ones measured through scanning electron microscopy (SEM) by Leung and Shapiro [26].

131
In the current work, we adopt Eq. (5) as a first approximation to represent the matrix wedge 132 which is to be spalled off later on. Then trial runs are performed to determine the moment at which the 133 matrix wedge is deactivated. From then on, the elements within the matrix wedge cease to contribute 134 to the overall stiffness. It needs to be pointed out that, such a treatment is only to save computational 135 time without the detriment to the numerical results, since the first principal stresses within the matrix 136 is checked and the spalled length is adjusted if necessary, to ensure that the matrix would not present 137 hyper-strength.
in which, d f is the fibre diameter, q is the inclination angle, P max is the peak pullout load of the aligned 129 fibre, whereas f t is the tensile strength of the matrix. The estimated values from Eq. (5) closely match 130 the ones measured through scanning electron microscopy (SEM) by Leung and Shapiro [26].

131
In the current work, we adopt Eq. (5)    by the matrix at the cracked surface, the matrix will spall. An estimated length of the spalled matrix can be obtained according to the work of Laranjeira et al. [34] as follows: in which, d f is the fibre diameter, q is the inclination angle, P max is the peak pullout load of the aligned 129 fibre, whereas f t is the tensile strength of the matrix. The estimated values from Eq. (5) closely match 130 the ones measured through scanning electron microscopy (SEM) by Leung and Shapiro [26].

131
In the current work, we adopt Eq. (5)   As mentioned before, we aim to explicitly model the matrix, the fibre and the interface in between. In a 2D setting, as shown in Fig. 2, both the fibre and the matrix are represented with four-node iso-parametric elements, whereas the interface is discretised into pairs of two-node line segments. The constitutive behaviour of the interface is governed by a cohesive law, as shown in Fig. 2b.
can be obtained according to the work of Laranjeira et al. [34] as follows: As mentioned before, we aim to explicitly model the matrix, the fib between. In a 2D setting, as shown in Fig. 2, both the fibre and the matr four-node iso-parametric elements, whereas the interface is discretised into segments. The constitutive behaviour of the interface is governed by a co Fig. 2b.
fibre, whereas f t is the tensile strengt  As mentioned before, we aim between. In a 2D setting, as shown four-node iso-parametric elements, segments. The constitutive behavio Fig. 2b. The above configuration is implemented in the commercial software ANSYS 12.0, using its APDL programming language. Regarding the constitutive laws, linear elasticity for the matrix and bilinear plasticity with isotropic hardening for steel fiber are adopted, respectively. Leung and Shapiro, 1999 In order to assess the effect of fiber yield strength on reinforcement efficiency in terms of maximum crack bridging force and total energy absorption, Leung and Shapiro [26] carried out pullout tests for steel fibers of different yield strengths and inclination angles at 0 • , 30 • or 60 • . All of the fibers are 0.5 mm in diameter and 22 mm in length. The pullout specimens are blocks of 25.4 mm× 12.7 mm × 9.5 mm in size, with an effective fiber length (equal to the embedment length in this case), L f , of 10 mm. The material parameters for the matrix, the fiber and the interface are given in Table 1, whereas the yield and tensile strengths of the four types of fibers are listed in Table 2. Additionally listed in Table 2 is the critical fiber length, i.e., the maximum embedded length for a fiber to be pulled out from a matrix without rupture [24]. It is related to the maximum shear stress, τ max , fiber tensile strength, f r , fiber cross-section, A f , and perimeter, p f , as follows

Experimental Setup of
It needs to be pointed out that this estimation is for aligned fibers only; in the case of inclined ones, this length is considerably smaller due to the contribution of fiber bending. Table 1. Material parameters for the matrix and the fiber given in [26].  Table 2. Yield and tensile strength of the four types of fibers tested in [26]; the corresponding L c is listed for a diameter of 0.5 mm, where (1) and (2) are calculated using the tensile strength, f r , and the yield strength, f y , respectively.

Identification of the Fiber-Matrix Interface Properties
In order to extract the interface properties from experimental pullout curves, we first introduce the concept of apparent interfacial shear stress, which is defined as: where τ(s) is the actual bond stress distribution at the interface and L e is the initial embedment length, the value of which is equal to L f . Assuming the fiber-matrix interface property is uniform, the maximum pullout load and peak frictional load are respectively calculated as: The values for P max and P f are determined from the pullout response of aligned fibers, as shown in Figure 3. Note that the values for P max and P f are averaged for the four types of fibers listed in Table 2 to obtain those of τ max and τ f c , as well as their standard deviations in Table 3. The critical slip for interfacial bond, s b c , 0.3 mm, is determined through trial and error, so that the first decaying branch of the numerical pullout responses, as demonstrated in Figure 3, should fall within the experimental range. As regards the friction coefficient given in Table 3, it is estimated according to the experimental results of Chanvillard [22], which was also adopted by Laranjeira et al. [27].  [26] and its bi-linear numerical approximation (red dotted line). P max and P f are the peak and the transitional pullout forces. Table 3. Extracted parameters of the fiber-matrix interface from the experimental data of Leung and Shapiro [26].

Numerical Model
The in-plane dimensions and boundary conditions to simulate the pullout tests performed by Leung and Shapiro [26] are illustrated in Figure 3. Note that within a two-dimensional plain stress framework, the fiber thickness, T f , is calculated through Equation (9) so that the contact area at the interface is the same as that of the original one. In the same way, the fiber height, H f , is determined via Equation (10) Additional described in Figure 4 are the boundary conditions imposed. Vertical displacements are prevented on the top and bottom sides, whereas horizontal movements are impeded on the left. The right end of the fiber is fixed in the vertical direction so that only horizontal movement is permitted. Note that the same boundary conditions are imposed for inclined fibers. The pulling process is carried out with intervals of 0.001 mm in the horizontal direction until 0.3 mm, followed with increments of 0.01 mm until the end.

Mesh Description
A typical mesh and detailed element distribution around the fiber for the inclination angle of 30 • are demonstrated in Figure 5. Note that the right end of the fiber leans on a matrix wedge, which will spall later on. For this particular case, the matrix and the fiber consist of 2198 and 154 four-node solid elements, respectively, whereas 289 two-node contact pairs are placed at the interface. The mesh sensitivity analysis performed to achieve a balance between the computational efficiency and accuracy is going to be presented in Section 4. 10 worth noting that this numerical model is mesh sensitive. To specify, contact is a highly nonlinear problem so the mesh density will influence the efficiency of calculation. Moreover, a coarse mesh will lead to higher stiffness of elements thus causing larger pulling forces. After trial and error (discussed later), this mesh is the overall best one to compromise between the computational efficiency and accuracy.

Numerical Results and Discussion
In this section, mesh sensitivity analysis is first conducted along the fiber transverse and longitudinal directions in order to determine the particular mesh to be employed for further studies. Second, the entire pullout load vs. slip displacement curves are extracted to compare with those obtained experimentally by Leung and Shapiro [26]. Third, the von Mises stress and the first principal stress evolutions are explored both for the fiber and the matrix. Finally, the pullout work is obtained.

Mesh Sensitivity Analysis
The mesh sensitivity analysis is carried out for the inclination angle of 30 • and the yield strength of 635 MPa (Type 3 in Table 2). Two kinds of mesh sensitivities are studied: the refinement in the transverse direction and along the axial direction of the fiber. The former is to check the capacity of the mesh in the fiber to bear bending moments, whereas the latter is to assess if the discretisation is fine enough to resolve the slip length. In the transverse direction, the fiber is split into one, two or four divisions (see Figure 6a); the corresponding load-displacement curves are plotted in Figure 6b. On the one hand, convergence is achieved as the mesh is refined. On the other hand, slight differences among the three curves are perceivable. This means that the mesh sensitivity in the transversal direction is rather modest. Meshes of two divisions across the transverse direction are employed for further studies.
Along the longitudinal direction of the fiber, four different element sizes are considered, 0.303 mm, 0.222 mm, 0.135 mm and 0.068 mm, which lead to 33, 45, 74 and 148 divisions along the 10-mm length (see Figure 7a); the corresponding pullout load vs. slip end displacement curves are depicted in Figure 7b. Note that as the mesh gets finer, the response is smoother, and convergence is achieved. Specifically, the same peak load is obtained for the two finer meshes. In order to keep a balance between the computational efficiency and accuracy of the results sought, the mesh size of 0.135 mm is selected for further studies. It needs to be emphasized that the mesh sensitivity in the longitudinal direction is more pronounced than that in the transverse direction. In addition, from Figure 7b, it is observed that the maximum pullout load was achieved at a slip displacement of 0.07 mm, which agrees with the statement of Morton and Groves [30], who claimed that this value should be of the order of, but less than half, a fiber diameter.

Validation against Experimental Pullout Load vs. Displacement Response
In order to verify the previously developed methodology, we compare the entire pullout curves with their experimental counterparts given by Leung and Shapiro [26]. This comparison is displayed in Figure 8 for fibers inclined at 30 • and 60 • with four different yield strengths given in Table 2. Note that both the peak loads and the general tendency are well captured; the numerical curves fall within the experimental range; in particular, the rising tail at the end of each pullout process is also reproduced.  [26] comparison: complete pullout curves for the four yield strengths given in Table 2; the fiber is inclined either at 30 • (left column) or 60 • (right column).
It needs to be pointed out that as the fiber yield strength increases, the numerical pullout curves exhibit more oscillations. This means that the relative slip between fiber and matrix is more unstable; thus, finer meshes are needed to reduce the noise observed in those curves for fibers of higher strength.

Stress Evolution within the Fiber
Taking Type 2 fiber (yield strength 469 MPa) with an inclination angle of 30 • as an example, the von Mises stress evolutions for several characteristic points within the fiber are examined. These points are the pullout end A, the embedded end D and two intermediate ones B (location of the matrix spalling) and C, as depicted in Figure 9. Note that at Point A, the first peak stress was obtained when the pullout load reached its maximum due to interface debonding. Then, after a slight decrease, this stress increased again until yielding at the slip end displacement of 3 mm. Similar peaks are observed for B and C at slip displacement of 0.3 mm and a second peak upon yielding at 2 mm for Point B and 5 mm for Point C, respectively. The second peak is attributed to the stress concentration due to the cusp formed by matrix spalling. This is the snubbing effect described by Li et al. [24]. In Figures 10 and 11, the first principal stress distributions in the fiber at different loading stages are plotted for Type 2 fiber inclined at 30 • and Type 3 fiber inclined at 60 • , respectively. Note that during the pullout process, there are stress gradients both in the transversal direction and along the longitudinal one, which indicates a clear bending contribution. Stress concentration is also observed at the fiber exit point. In addition, the maximum tensile stress is always inferior to the fiber tensile strength. This means that the fiber was pulled out, but not broken. Furthermore, increasing separation is observed at the upper interface as the pullout process advances; this requires a numerical model that is able to capture a mixed-mode interface failure to correctly reproduce the experimentally-observed phenomenon.
There is an interesting result that goes against our intuition, which is that the stress along the fiber should be compressive on the concave surface and tensile on the convex one. However, the opposite is observed particularly in Figure 11. This is because a fiber segment bends the most right before passing the cusp; then, it is straightened by the pulling force, leading to compression on the convex part and tension on the concave part. Additionally observed from Figure 11 is that the length of the red region in the fiber is almost constant since the segment on the right gets unloaded as the left part is being stressed. This further demonstrates the stable deterioration of the interface during pullout. Only when a small fragment of the fiber is remaining in the matrix block, it serves as a hook, which grips the matrix tightly due to its relatively large stiffness. This leads to the tail-rise phenomenon, which is increasingly obvious when the yield strength is larger, since a stronger grip is formed.    Figure 11. The first principal stress evolution for pullout displacement from 0.3 mm to 7.8 mm for a fiber inclination of 60 • and a yield strength of 635 MPa (Type 3 in Table 2).

Stress Evolution in the Matrix
For the matrix, we are more concerned with the tensile stress distribution to ensure that no fracture should take place where matrix spalling is not expected. Three representative points, E, F and G (see Figure 12), are selected to display the first principal stress evolution in the matrix. Point E is where the matrix is expected to spall. Point G is the location where the fiber is anchored, whereas F is the point in the matrix close to the fiber center. Since the tensile strength was not measured, we estimate it to be 1/12 of the compressive strength, which is 3.0 MPa, as shown in Table 1. Note that at Point E, there is relatively large stress oscillation during the pullout process. This indicates, on the one hand, the stressing-relaxing cycle endured by the matrix. On the other hand, it can be attributed to the fact that the spalled matrix is assigned with a zero stiffness at the moment of spalling, whereas the real failure process is gradual. The stress evolution curves at Points F and G, assimilate those of global pullout curves in Figure 8, each with a different amplitude. Furthermore, it is noted from Figure 12 that the maximum tensile stress due to the axial pullout of the fiber and bending load occurs at the close region of the fiber exit point. This agrees with the assumption adopted by Zhang and Li [37] in their study on the effect of inclination angle on fiber rupture load in fiber-reinforced cementitious composites.

Variation of the Pullout Work with Respect to Fiber Yield Strength
The pullout work is calculated as the area under the load vs. slip displacement curve. In Figure 13, both experimental and numerical values for fibers inclined at 30 • and 60 • are depicted. Note that as a general trend, the pullout work increases with the increase of fiber yield strength, and such a tendency is correctly captured by our numerical model.

Parametric Studies
In this section, the influences of both material properties (fiber yield strength, matrix tensile strength, interface bond) and geometric factors (fiber diameter, embedment length and inclination angle) on the overall pullout behavior, in particular on the maximum pullout load, are studied.

Interface-Related Parameters
The influence of the bond strength, the critical slip displacement, the internal frictional resistance and the friction coefficient is herein explored.

Interfacial Bond Strength and Internal Frictional Resistance
As described above, the interfacial bond between fiber and matrix is mainly composed of two parts, internal frictional resistance and bond strength. Simulations focusing on varying τ max and τ f c separately for fibers inclined at 0 • , 15 • , 30 • and 60 • are carried out. The obtained pullout curves are demonstrated in Figure 14, where three different values, 1.3, 2.7 and 5.4 MPa for τ max , 0.25, 0.50 and 1.0 MPa for τ f c , are considered. It is noted that τ max has a dominant influence at the pre-peak stage (when the slip-end displacement is less than 0.3 mm), in particular over the peak load, whereas τ f c plays a more important role during the post-peak phase. In addition, the total slip displacement decreases as the inclination angle increases. This means the critical embedment length decreases with the inclination angle.
It needs to be pointed out that the highest peak load is obtained when the inclination angle is at 15 • for all cases. This is confirmed in Section 5.2.5, in Table 5, where θ max falls between 12 • and 15 • for different values of fiber yield strength.

Critical Slip for Interfacial Bond
As mentioned before, the critical slip displacement, s b c , denotes the maximum value for which interfacial bonding due to internal physical and chemical cohesion is active. Three different values, 0.15, 0.30 and 0.60 mm, of s b c are selected to demonstrate such an effect for fibers inclined at 0 • , 15 • , 30 • and 60 • , respectively, in Figure 15. Note that the difference of the pullout curves is only perceived at the first decaying branch; larger s b c leads to a less steep drop from the maximum pullout load.

Friction Coefficient
As mentioned before, the interfacial friction consists of an internal friction and a pressure-dominated one due to fiber inclination (dry friction). The former is considered as a property of the interface, whereas the effect of the latter is only revealed when the fiber is inclined to the load direction. Four inclination angles (0 • , 15 • , 30 • and 60 • ) and three different values of the friction coefficient (0.3, 0.6 and 0.9) are simulated. The obtained pullout curves are demonstrated in Figure 16. Note that a higher friction coefficient leads to a load increment during the entire pullout process, and this effect is more pronounced for larger slip displacements. This evolution is attributable to the snubbing effect around the matrix cusp and corresponds well with the results of Cailleux [29].

Fiber-Related Parameters
Parameters related to fiber, such as the yield strength, fiber diameter and length and the inclination angle, all play an important role in the pullout curves. Herein, the influence of each parameter is studied while keeping the remaining ones unchanged.

Fiber Yield Strength
The influence of fiber yield strength has been experimentally demonstrated by Leung and Shapiro [26] and our numerical results in Figure 8. It is expected that when a fiber of the same size and surface treatment is pulled out, the same amount of energy will be consumed to overcome the interface bonding and frictional resistance. Furthermore, since the coefficient of friction is mainly dependent on the surface property, it is assumed constant for the same type of matrix and fiber. Consequently, the frictional energy will vary little with the fiber yielding strength. However, more energy is needed to bend the fiber as the yielding stress increases. This is clearly illustrated in Figure 17, where the entire pullout curves are moved upwards as the fiber yield strength is augmented. An interesting observation from Figure 17 is that the tail of the curve goes up, which was also seen in the experimental curves of Leung and Shapiro [26]. When the fiber is nearly pulled out, the short remaining end is rather stiff; thus, it may rotate as a short rigid cylinder, which subsequently cuts into the matrix. This has been explained previously through the first principal stress distribution in Figure 11, where the red areas represent the regions in heavy tension (nearly yielding), whereas the blue ones stand for the regions under less tension or compression.

Fiber Diameter
The overall pullout response is closely related to the fiber diameter. First, a larger diameter will lead to a wider contact area between the fiber and the matrix. Consequently, a higher pullout load is attained. Second, different diameters lead to different moments of inertia of the fiber cross-section; thus, the dissipated energy will differ. To explore this influence in detail, three values of fiber diameter (0.25 mm, 0.50 mm, 1.00 mm) with four inclination angles (0 • , 15 • , 30 • and 60 • ) are selected to perform the numerical pullout simulations. The related parameters are documented in Table 4, whereas the pullout responses are shown in Figure 18. Additionally given in Table 4 is the aspect ratio, L e /d f , which is herein defined as the ratio between the fiber embedment length, L e , and the fiber diameter, d f . It is commonly accepted that the higher the aspect ratio is, the more efficient is the fiber at bridging cracks.  It is observed from Figure 18 that it is more difficult to pull a thicker fiber out. Consequently, more energy is consumed to pull an inclined thick fiber. Furthermore, the phenomenon of the rising tail of the pullout response is particularly noticeable for thicker fibers. In addition, the effect of bending is more pronounced for a thicker fiber, whereas more energy is consumed at the interface for a slim fiber. This explains the result that a larger peak pullout value is obtained at 30 • than at 15 • for the fiber of 1.0 mm in diameter, whereas the opposite is observed for the fibers of 0.25 mm and 0.5 mm in diameter. At the same time, it is noted that when a thin fiber is pulled, it may yield or get ruptured instead of being pulled out. This is the reason that the curve fluctuates around 5-mm-slip in the case of the 0.25 mm-diameter fiber inclined at 60 • .

Fiber Length
In Figure 19, the pullout responses for fibers of different lengths from 4 mm to 17.5 mm inclined at 0 • and 30 • are illustrated. Since a longer fiber length also means a wider interface area, thus a higher pullout load is attained. For a more clear representation, the apparent interfacial shear stress vs. relative slip normalized by the initial embedment length, which is equal to L f , is also shown. It is observed that, for both aligned and inclined fibers, convergent results are achieved for longer fibers. This indicates that, when the rest is kept the same, the apparent bond strength coincides with the real one if the fiber embedment length is relatively large. Otherwise, the interface strength obtained through the pullout tests will be over-estimated.

Fiber Inclination Angle
The effect of the fiber inclination angle on the entire pullout response is explored for six orientations from 10 • to 70 • as demonstrated in Figure 20. Note that peak values for fibers inclined at 10 • and 20 • are higher than the rest, whereas the post-peak values during 2 to 7 mm only vary slightly. This is due to the fact that the pullout force is only a component of the bonding force before de-bonding. Afterwards, the contribution of fiber bending and friction is increased for all inclination angles. 26 Furthermore, several inclinations (10º, 20º, 30º, 45º, 60º, 70º) are simulated to explore the influence of spalling length on the whole curve of pullout load-displacement.  Figure 21, it could be seen that less inclination leads to smaller peak load. This is because the pullout force is just a component of the bonding force before debonding.
Afterwards, the contribution of fibre bending and friction turns greater for all inclinations. The post-peak values during 2-7mm are almost insignificant for all inclinations within the variation of 10N while the fibre of 10º is the lowest. From 10º to 45º, with the increment of angles, the post-peak value climbs up as the ones for larger inclinations of 45º, 60º, 70º stays nearly constant. This means the greater bending caused by larger inclination is compromised by matrix spalling.

The effect of spalled length
As discussed previously, matrix spalling length is vital while assessing the crack-bridging capability of fibres because it will significantly affect the bending

The Maximum Pullout Load
To explore the effect of fiber inclination, a spectrum of angles up to 85 • is simulated by keeping the fiber, matrix and interface properties constant. It is known that the length of the spalled matrix and the moment when the matrix spalls both matter in the pullout responses. According to Laranjeira et al. [35], spalling is considered to take place right after the beginning of fiber debonding, but prior to its full accomplishment. After some trial runs, this slip displacement is estimated, which is around 0.01 mm. The spalling length is calculated according to Laranjeira [35], which corresponds well with that of Leung and Shapiro [26]. Simulations of different fiber yield strengths are carried out, and the obtained pullout curves are plotted in Figure 21. The fitted function for the maximum pullout load in terms of fiber inclination angle is given by Equation (11), and the corresponding fitting parameters are listed in Table 5. Table 5. Fitted parameters for the maximum pullout load vs. fiber inclination angle, where θ max is the angle where the maximum peak load is attained.   Table 2 and fitted curves using Equation (11) with parameters given in Table 5.

Fiber Type
From Figure 21, when the fiber is inclined at θ max , which is around 15 • , the maximum pullout load is the largest. After that, the maximum pullout load goes down almost linearly. This differs from the result of Morton and Groves [30] for a polyester resin matrix. This indicates that the optimum inclination angle for maximum pullout resistance varies with both fiber and matrix strength, as well as the interface properties.

Matrix-Related Parameter: The Tensile Strength
The influence of matrix strength on the fiber pullout responses has been actively studied in the last three decades [24,[38][39][40]. However, since the change in matrix strength may inevitably lead to slight [24] or marked [38][39][40] variation in the interface properties, the experimental observations can differ significantly. In the current numerical study, we isolate the matrix strength from the rest of material and geometric parameters to study its effect on the global pullout response. This is carried out through its direct influence on the matrix spalling length [35]. Taking the fiber with a yield strength of 635 MPa as an example, pullout out responses for three levels of matrix tensile strength (1.5, 3.0 and 6.0 MPa) and fiber inclined at 0 • , 15 • , 30 • and 60 • are explored. The corresponding spalling lengths calculated according to the work of Laranjeira et al. [35] are listed in Table 6. The obtained pullout curves are plotted in Figure 22. Note that a stronger matrix leads to a shorter spalling length and results in a higher peak load and a larger amount of absorption energy. However, a small difference is perceived for slip displacements between 3 mm and 7.5 mm for the fiber inclined at 30 • , whereas more variations are observed for fibers inclined at 60 • .

Conclusions
We have proposed a numerical model to explicitly reproduce the pullout behavior of a single fiber embedded within a cement-based matrix. This model takes into consideration the gradual deterioration of the interface bond, internal and dry friction, as well as matrix spalling. In particular, a constitutive law, which isolates the contributions of the internal bond, internal friction and dry friction is formed and validated. Cohesive elements endorsed with mixed-mode fracture capacities are implemented to represent the bond-slip behavior at the interface. Contact elements with Coulomb's friction are placed at the interface to simulate frictional contact. Matrix spalling is modeled through material erosion to save computational time. The bond-slip behavior is first calibrated through pull-out curves for fibers aligned with the loading direction, then validated against experimental results carried out by Leung and Shapiro for steel fibers oriented at 30 • and 60 • . The influence of fiber yield strength on the stress distribution within the fiber and the matrix, the effect of the inclination angle, as well as the matrix spalling length on the pullout response are all explored in detail.
The proposed methodology provides the necessary pull-out curves for a straight fiber oriented at a given angle for multi-scale models to study fracture in fiber-reinforced cementitious materials. The novelty lies in its capacity to capture the entire pullout process for a fiber with an arbitrary inclination angle. It is equally applicable for hooked or deformed fibers with increased computational efforts.
In addition, the following specific conclusions can be drawn from the parametric studies: • Fiber critical length diminishes considerably for larger inclination angles.
• The matrix strength has a dominant influence at the initial stage of pullout, particularly on the maximum pullout load. • The inclination angle at which the maximum peak pullout load is obtained differs when fiber, matrix and interface properties are different. Author Contributions: Both authors closely interacted to conceive the idea and the structure of the paper. Hui Zhang carried out the numerical simulations and provided the first draft. Rena C. Yu supervised the work and re-edited the final version.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: