A micromechanical model of fiber bridging including effects of large deflections of the bridging fibers

Abstract A micromechanical model of cross-over fiber bridging is developed for the prediction of macroscopic mixed-mode bridging laws (traction-separation laws). The model is based on non-linear beam theory and takes into account debonding between fiber and matrix as well as buckling of fibers in compression. Further, it is shown how failure of the bridging fibers can be taken into account through a Weibull distributed failure strain. Predictions made by the proposed model are compared with predictions made by detailed 3D finite element models, and a very good agreement was observed. It is shown that models based on linear beam theory are only valid for small transverse deflections of the bridging ligament and greatly underestimate the force transferred by ligaments subjected to moderately large deflections. The novel model, on the other hand, is applicable in the entire range where the bridging problem transitions from a beam bending problem to a bar-like problem. Finally, an example of how the proposed model can be used for parameter/sensitivity studies is given. A conclusion from this study is that reducing the fracture toughness, G c , of the interface between fibers and matrix may lead to increased energy dissipation through cross-over fiber bridging as more fibres remain intact longer.


Introduction
In some fiber reinforced polymer laminates, the fracture process zone (FPZ) of a delamination can be long in comparison with laminate dimensions due to the development of fiber bridging in the wake of the crack tip as illustrated in Fig. 1. Fibers that bridge the fracture surfaces transfer tractions between the two surfaces (see Fig. 2) and can enhance the delamination resistance substantially as the crack extends [1][2][3][4][5][6]. Micromechanical models of cross-over fiber bridging can be valuable tools for studying the underlying mechanisms and how to utilize fiber bridging to maximize the fracture resistance and therefore the damage tolerance.
A number of micromechanical models for the prediction of macroscopic traction-separation laws for cross-over fiber bridging have been developed. Spearing and Evans [1] developed a model of cross-over fiber bridging in pure mode I delamination including shear deformations in a bridging ligament with rectangular cross section [7]. Shear deformations will dominate in short ligaments, i.e. at small opening displacements. At large opening displacements the bridging ligament becomes slender as it peels off the fracture surface, and shear deforma-tions become negligible. Kaute et al. [8] proposed a model where only the axial stiffness of the bridging fiber is considered. They assumed the fiber to be straight and long in comparison with its diameter. A reduction in tensile strain in the fiber due to slipping of the fiber within the uncracked matrix was taken into account through fracture-mechanics considerations. The model predicts that the normal tractions acting on the fracture surface from a single fiber will increase to a plateau for increasing opening displacements. A length-dependence of fiber strength was included through a Weibull distribution. In this way, the number of surviving fibers decreases for increasing opening displacements and the resulting traction on the fracture surface also decreases, as observed experimentally. Ivens et al. [9] developed a model based on the work of Wells [10]. A simple model of a DCB specimen with one fiber bridging the crack was analysed on the basis of energy considerations. The fiber was assumed to be straight and only transfer longitudinal forces. More recently, Daneshjoo et al. [11] published a model building on the work by Kaute [8]. Their model considers matrix spalling, fiber pull-out and fiber fracture as the main failure mechanisms. The absorbed energy in the fiber bridging zone is obtained by summing the energy terms associated with each of these mechanisms. Sørensen et al. [6] proposed a micromechanical model for mixed mode delamination based on classical Euler-Bernoulli beam theory. Their work can be viewed as an extension of the model by Spearing and Evans to mixed mode I/II delamination. The two models are identical for pure mode I opening if the shear term of Spearing and Evans is omitted. For that case both models predict the normal traction to be inversely proportional to the square root of the normal opening. Both these models are limited to infinitesimally small deflections of the bridging ligament, i.e. when the local normal opening displacement is much smaller than the height of the bridging ligament. The models proposed by Kaute et al. [8], Ivens et al. [9] and Daneshjoo et al. [11], on the other hand, are only applicable when the local normal opening displacement is several orders of magnitude larger than the height of the bridging ligament. A comprehensive review of fiber bridging investigations, including micromechanical models for crossover fiber bridging in the wake of a delamination front under different modes, is given by Khan [12].
The scope of the present work is to establish a micromechanical model applicable to the full range of deformations that the bridging ligaments are subjected to in mixed mode I/II delamination. I.e. the model should be able to predict the transition from a beam bending problem to a bar-like problem. The model is based on the framework of Sørensen et al. [6], extended with moderately large deflection beam theory (including von Kármán strains) and a Weibull distributed failure strain of the fibers. Fig. 3 shows schematically a fracture process zone with a single bridging fiber. This corresponds to an arbitrary position within the fracture process zone depicted in Fig. 1. The fiber is considered as a beam. In the depicted case, the beam is fixed at the left end while the local opening displacements are imposed to the right end. The horizontal and vertical displacement at the right end are denoted δ x and δ y , respectively. These correspond to the local opening displacements which are the prescribed input to the model. Both ends are constrained from rotating. The prescribed end-displacements result in varying displacements along the beam. The horizontal displacement as function of position is denoted uðxÞ and the corresponding vertical displacement is denoted wðxÞ. The solution of this beam problem forms the basis for the micromechanical model that will be described in the following subsections. In short, the tractions acting on the fracture surface from one single ligament are calculated from moderately large deflection beam theory. Buckling is considered for fibers subjected to compression, and limits the contribution from compressed fibers. The details of the beam model are described in Section 2.1. The length of the bridging beam is determined from Griffith-like energy considerations, as described in Section 2.2. Finally, the statistical contribution from a large number of bridging fibers is included in Section 2.3 where also fibre failure is taken into account through a Weibull distributed failure strain.

Micromechanical model
A Matlab implementation of the model has been made and results obtained with it will be compared to finite element predictions in Section 3.

Determining forces from a single fiber with given length
We assume finite rotations and deflections in the following. The material of the fibers is assumed to be linear elastic, and the laws of  Fracture process zone with bridging tractions. Both the normal traction, σ n , and the tangetial traction, σ t , are functions of the normal opening, δ n , and the tangential opening, δ t . Fig. 3. Idealized fracture process zone with a single bridging fiber. elasticity remain the same as for classical beam theory. Thus, the axial force N and bending moment M can be written where E is the Young's modulus, A is the cross sectional area, I is the second moment of area, ε n is the strain at the neutral axis and κ is the curvature. In the case where the bridging ligament consists of a single fiber, the elastic modulus E equals the fiber modulus E f . In case the beam represents a ligament consisting of several fibers and matrix, the modulus is a representative Young's modulus of the composite material, E c . In comparison with classical beam theory, the expression for the axial strain at the neutral axis contains a second term due to finite rotations (von Kármán strain [13]) where u is the horizontal displacement and w is the vertical displacement. The definition of curvature also has a nonlinear rotation term [14] κ ¼ À However, we will restrict our model to moderately large rotations where the square of the slope is small compared to unity and the curvature can be defined in the same way as in the theory of small deflections Equilibrium in the vertical direction for a beam with axial force, but no distributed vertical load, is given by [15] À EI where both the axial force, N, and the vertical displacement along the beam, wðxÞ, are unknowns that must be determined. We assume the axial force as well as the cross-sectional properties to be constant along the beam. By dividing Eq. (6) by EI and introducing λ 2 ¼ N EI one gets Depending on the axial force, the general solution of this ordinary differential equation is [16] wðxÞ The constants C 0 to C 3 as well as λ (i.e. five unknowns) must be determined from the boundary conditions. There are six boundary conditions in Eq. (9), but both uð0Þ and uðLÞ are required to determine ΔL and thus λ. Note that the constants will differ for the two cases (N ¼ 0 and N -0). The boundary conditions of the beam shown in Fig. 3 are If the axial force is kept constant for now, then Eqs. (8) and (9) (the parts related to w) form a system of four linear equations that must be solved simultaneously and the following values can be determined for the four constants Since the axial force is assumed to be constant along the beam, it can be expressed in terms of the total elongation of the beam.
As can be seen from Eq. (3), the axial strain depends on both vertical and horizontal displacements, u and w. The elongation of the beam can be determined from the integral of Eq. (3) which equals the difference in horizontal movement at the ends plus the integral of the axial strain araising from vertical movement along the beam Inserting (12) into (11) yields Inserting (8) into (13) results in rather complicated transcendental equations for λ, from which no closed form solution has been obtained. Instead, the axial force was determined iteratively. The axial force in iteration (i + 1), N iþ1 , was found by inserting the axial force from iteration i, N i , into Eqs. (8) and (10) and using this in Eq. (13). The vertical displacement wðxÞ calculated from linear theory (i.e. setting N ¼ 0) can be used as a starting point for calculating the axial force iteratively. Combining Eqs. (8) and (10) gives the solution for N ¼ 0 Differentiating Eq. (14), inserting into Eq. (13) and performing the integral gives This initial estimate is normally very good and only few iterations are needed to obtain acceptable accuracy.
Since the ends of the beam are constrained from rotating, the tangential force acting on the fracture surface from one beam is identical to the axial force N. Similarly, the normal force is identical to the shear force V at the end of the beam Slender fibers are susceptible to buckling if they are subjected to compressive axial forces approaching the critical load, which for fixed ends can be written [15] N cr ¼ À Both the tangential and the normal force are therefore set equal to zero if N < N cr .

Determining the undeformed length of the bridging fiber
The forces in a bridging beam with a given undeformed length L were derived in the previous subsection. The length of the beam has to be determined for given opening displacements. Assuming that the fracture process zone of the debonding process between fiber and matrix is small, the length can be determined by requiring that the energy release rate of the bridging mechanism equals the fracture energy of the interface G c . The potential energy of an elastic body, Π, is defined as follows [17] where U is the strain energy stored in the body and F is the potential of the external forces. For the present problem we prescribe displacements not forces, so the last term in (18) vanishes, so F ¼ 0 and Π ¼ U. Therefore, the rate of change in potential energy with the crack area, G, can be written where A is the fracture area and b is the effective width of the fracture area.
Since the material is assumed to be linearly elastic, the strain energy density can be calculated aŝ where σ and ε denote stress and strain, respectively. Then the total strain energy can be expressed as where y is the distance to the neutral axis. The first term is due to stretching and the second due to bending. The bending moment M is given by MðxÞ ¼ ÀEIw 00 ðxÞ where wðxÞ is the deflection of the nonlinear beam as described in the previous subsection. Numerical differentiation of the strain energy with respect to the length of the beam and Ridders' method is used in our Matlab implementation to determine L so that G ¼ G c by solving Fig. 4 shows two fibers bridging the fracture surfaces in opposite directions. As can be seen, the contribution from both the normal and shear force to the fracture surface tractions is depending on fiber orientation.

Contribution from a large set of fibers
Let η u and η d be the number of fibers bridging a unit area of the fracture surface in the upwards and downwards direction, respectively, see Fig. 4. The resulting tractions on the fracture surface can then be written where σ t and σ n are the tangential and normal tractions, respectively. The sign convention follows from Fig. 4.
If the number of fibers initially bridging the fracture zone is very large, then the fraction of active fibers can be approximated by a continuous function where η u0 and η u0 are the amount of fibers initially bridging a unit area in the upwards and downwards directions, respectively. The fraction of active fibers is represented by f u and f d for the two diagonal directions. It is necessary to keep track of the number of fibers crossing each direction, since it can be different. A mode II opening component will cause tension in one direction and compression in the other. These stresses/ strains may be superimposed to the tension resulting from a mode I opening component. Therefore, the stress/strain level may be different in fibers bridging in opposite directions. Note that Eqs. (23) and (24) relate the forces in the representative volume element (RVE) at the micro scale to the statistically homogeneous tractions at the macro scale. This is a critical step that is denoted homogenization in micromechanics [18].
For brittle fibers, such as carbon or glass, the strength is normally limited by the most severe defect present. Literature suggests that the strength of individual fibers and fiber bundles subjected to uniform tension are well described by the Weibull distribution [19][20][21]. Under this assumption, the fraction of still intact fibers can be expressed as unity minus the cumulative distribution function of the Weibull distribution [21] f ¼ e where ε is the applied strain and ε 0 and m represent the Weibull scale and shape parameters for failure strain. ε 0 can be thought of as the strain associated with a probability of failure of 0.63 for a length L 0 of fiber strained uniformly, while m describes the flaw distribution and the size dependency. Except for the case of pure mode II opening displacement, the fibers will not be subjected to a uniform strain. Instead, the strain will vary both along and across the fiber. Assuming that fiber failure is initiated at surface defects, Sørensen and Goutianos [22] proposed a surface integral to be used for fibers subjected to a heterogeneous strain field where εðx; θÞ is the current strain at the fiber surface at position ðx; θÞ and D is the diameter of the fiber. It should be noted that the Weibull distribution is only defined for ε ∈ 0; 1 ½ i. Since it is also believed that brittle materials are more prone to fail from tension than compression, we have chosen to only include the contribution of tensile strains to the probability of fracture and hi denotes Macaulay brackets in Eq. (26).
The largest combination of bending and axial strains occurs at the ends of the free span of the partially pulled off fiber (the anchoring point). In this way, the cross-sections at the ends account for the major contribution to the surface integral in Eq. (26). However, as the fiber continues to be pulled off, the end of the free span translates along the fiber and the previously most strained cross-sections are partially unloaded. This causes Eq. (26) to reach a plateau when the opening mode is constant. It therefore fails to capture size (length) effects whereby the probability of failure due to severe defects increases as an increasing length of the fiber has been subjected to severe loading. Ideally, one should consider the maximum strain that each point in the fiber has experienced. However, it is not possible to derive any analytical expression for the maximum strain a point has experienced since the shape of the fiber is determined iteratively for each opening state. Instead, it was chosen to integrate strains along the fiber's endcircumference in the current state and multiply by the peeled off length of the fiber as follows (as if the entire free-spanning fiber length has experienced the same strain as the current anchoring point): In pure mode II, the strain is constant in the entire fiber regardless of the magnitude of the opening displacement and Eq. (27) captures the size effect correctly. However, when there is a mode I component, the curvature at the end increases with increasing opening displacement and the current strain at the end will be greater than the strain at the end in previous opening states. It is still believed that Eq. (27) is more representative than Eq. (26).
The axial strain at any point in the fiber can be expressed as At the end of the fiber the axial strain reduces to

Validation of the model
The model described in the previous section was implemented in Matlab. Comparisons against finite element analyses (FEA) will be used to verify the implementation and assess the accuracy of the assumptions and simplifications made in the proposed micromechanical model. A set of representative, but somewhat arbitrary, material and geometry parameters was chosen to test the model. These parameters are summarized in Table 1.
Finite element analyses (FEA) of the bridging problem were carried out using the non-linear finite element code LS-DYNA. Fig. 5 shows the parts in the model as well as the mesh of the cross-section of the fiber. A structured mesh of approximately 810,000 brick elements (eightnode hexahedron) was used to represent the fiber (red in the figure). This very fine discretization was chosen in order to have a sufficient resolution of the fracture process zone between fiber and matrix. Eight-node cohesive elements were used to model the interface between the fiber and the matrix (blue in the figure, seen through the partly transparent matrix parts). A bi-linear cohesive law with an initial stiffness of 20 GPa/μm, a maximum traction of 4 GPa and a critical separation of 0.5 μm was used. The critical separation corresponds to D/20, so the size of the active fracture process zone is small and thus satisfies LEFM conditions, i.e. the dissipative processes remain confined to a region in the vicinity of the crack tip that is small compared to the structural dimensions. The exact same traction separation law was assumed for mode I, II and III. Mode interactions were modelled using the power law criterion [23] with the exponent equal to 2, i.e. quadratic interaction. This means that the direction of the traction vector follows that of the opening vector. Since we assume the same traction separation law for all modes, the combined work of the cohesive tractions doesn't depend on opening path history [24]. The matrix was modelled with two rigid bodies, representing the two fracture surfaces. The lower surface (green in the figure) was fixed while opening displacements were imposed on the upper surface (yellow in the figure).
A comparison of the normal force predicted by FEA, the novel nonlinear model as well as the model proposed by Sørensen et al. [6] (modified to circular cross-section) for a single fiber bridging upwards in pure mode I opening is shown in Fig. 6. The normal opening, δ n was increased incrementally for the two analytical models. The peeled off length, L, and the forces acting in the fiber were determined for each state as described in the previous section. Very good agreement between the non-linear model and FEA is observed, while the model based on linear beam theory underestimates the force level significantly for opening displacements exceeding the fiber diameter. Although the force comparison is a very good indication that both the free span, L, and the slopes and curvatures along the fiber, wðxÞ, are well predicted, additional comparisons are given in Figs. 7 and 8. The finite element model predicts the crack front between the fiber and the matrix to be curved. It is difficult to give a stringent definition of L in this case. Both the maximum and minimum value depending on where it is measured are given in Fig. 7. Regardless of which of these are chosen as the reference, the proposed model predicts the free span to be slightly longer. This has an impact on the predicted slope as seen in Fig. 8, which shows the deflections predicted at one of the last converged states of the FEA with the proposed model prediction overlaid. The proposed model predicts the rotation at the mid-span to be slightly smaller than what predicted by the finite element model, but the overall agreement is good. Note that at this stage, the non-linear model predicts the fiber to act more like a bar than a beam.
The mode-mixity can be defined as φ ¼ atanð δt δn Þ where δ t and δ n are displacement jumps between two points located on opposite crack faces at some distance behind the crack tip. Four mode mixities were selected for validation of the non-linear model against FEA, namely φ ¼ 0 ; φ ¼ 45 ; φ ¼ 67:5 and φ ¼ 90 . The results are shown in Fig. 9 where the total opening is defined as δ ¼  bridging both diagonal directions, and the two curves are therefore on top of each other in Fig. 9(a). Fig. 10 shows the bridging tractions predicted by the proposed micromechanical model for the chosen mode angles assuming η u0 ¼ η d0 ¼ 1. Except for the pure opening modes, the ratio of the normal to tangential traction will vary with the magnitude of the opening. Thus, the direction of the traction vector will vary and in general will be different from the direction of the separation vector. The mixed mode cohesive laws are thus coupled. These findings are in agreement with the finding of the simpler model [6].
Small rotations were assumed in Eq. (5). It is important to check the validity of this for the bridging problem. The value chosen for G c in Table 1 corresponds to 1000 J/m 2 . This value is significantly higher than what is typically reported in the literature. This was initially believed to be conservative for the purpose of validating the model. However, the predicted slopes of the fibres are smaller than the maximum seen in delamination tests. The fracture energy of the fibre/matrix interface was therefore varied to investigate the effect of this parameter on the predicted maximum slope of the fibre. The rotation at the mid-span is defined as: Fig . 11 shows the predicted rotation of the fibre at the mid-span for various mode I opening states for a wide range of G c . A very good agreement between FEA (circles) and the micromechanical model (continuous lines) is seen for the reference G c (1000 J/m2) and lower values of G c . The predictions for the quadruple G c (4000 J/m2) are more apart, but the agreement is still acceptable. It should be noted that values stated in published literature is in the range 10-300 J/ m 2 [25][26][27][28][29][30]. Still, rotations grater than 25°are frequently observed. The reason for this is likely a deeper embedment than half the diameter and that energy is dissipated in the bulk of the matrix material and not only at the fibre/matrix interface.

A short parameter study
The proposed semi-analytical model has been used in a limited parameter study to demonstrate its usefulness for such. The parameters varied were the fracture energy of the fiber/matrix interface, G c , and the fiber diameter, D. Note that in the case where D was scaled, so was η to reflect a constant fiber volume fraction. The values of the parameters used in this study are given in Table 2. The remaining input parameters are as given in Table 1.
For simplicity, only pure mode I and II delaminations were investigated in this limited study and buckling of compressed fibers in mode II was neglected. The predicted bridging laws are shown in Fig. 12. The energy dissipated through the bridging mechanism can be found by integrating the traction separation curve. The resulting curves are shown in Fig. 13. As can be seen, increasing G c leads to a more rapid   decrease in tractions and therefore the dissipated energy is reduced. This effect is stronger in mode-I than in mode-II.
It is possible to derive analytical expressions for the tractions and dissipated energy in pure mode-II delamination. The elastic energy in a fiber of linearly elastic material with diameter D and length L subjected to an elongation δ t can be written Differentiating with respect to L gives Solving for the peeled off length assuming that half the circumference was initially embedded in matrix: The strain in the fiber then becomes I.e. the strain (and stress) is constant regardless of opening displacement in pure mode II. The force in a single bridging fiber can be written The bridging tractions then become The fraction of intact fibers can be written where the relation εL ¼ δ t has been used and Φ has been introduced as: The energy dissipated per unit area can be found by integrating tangential tractions with respect to tangential displacements. The energy required to completely separate the two surfaces can then be found by integrating to infinite displacements From the expression above, it can be seen that for m ¼ 2; G c will cancel out. For m > 2; W will decrease for increasing G c . For m < 2; W will increase with increasing G c . It can also be seen that W / D m=2 . However, η 0 is likely to be a function of D. If the fiber volume fraction is kept constant, then the number of fibers, η 0 , will be inversely proportional to the square of the diameter, D, i.e. η 0 / D À2 . In this case, it can be argued that W / D m=2À2 . Then D will cancel out for m ¼ 4. For lower values of m; W will decrease with increasing D while it will increase when m > 4.

Discussion
The results presented in Section 3 clearly demonstrate the need to take geometric stiffness effects into account for accurate modelling of large-scale bridging problems. As can be seen in Fig. 6, the linear model becomes inaccurate for mode I when the deflection approaches a magnitude of the same order as the fiber diameter and greatly underpredicts the force level for large opening displacements. As can be seen in Fig. 9(a), the axial force is dominating the response of the fiber in pure mode I already at a normal opening displacement equivalent to half the fiber diameter. Therefore, a linear beam model is incapable of capturing the real physics of the problem already at small normal opening displacements. However, as the last term in Eq. (12) and (13) become zero in pure mode II, the proposed model and the linear model produce identical predictions in this case.
The finite element model includes large rotations and shear deformations. The very good agreement observed between the FE model and the proposed analytical model when it comes to force predictions indicates that (at least for the chosen set of parameters) the proposed model has sufficient accuracy even if it is constrained to moderately large rotations and the Kirchhoff hypothesis (plane cross-sections remain plane). Neglecting the square of the slope in the denominator of Eq. (4) introduces an error. The maximum rotation of the centre point of the fibre was 18.2°for the reference case (see Fig. 11). This would lead to an error of 16.7% for this particular point along the fibre. At the same time, the numerator is zero for this particular point. The fact that the maximum curvature occurs at the ends where the rotation is zero, while the curvature is zero at the mid span where the maximum rotation occurs, may explain why the model predicts seemingly accurate results even when the centre rotation is quite large.
Rotations are well within the range of moderately large deflection beam theory when G c is chosen in the range reported in the literature on testing adhesion between fibre and matrix [25][26][27][28][29][30]. However, larger rotations are observed in delamination tests. This can only be explained by fibres being over-embedded and mechanically licked by the matrix material.  The analytical model predicts a slightly longer peeled off length L than the FE model and thus that the rotation at the mid-span is slightly smaller. It should be noted that the boundary conditions (BCs) differ slightly between the FE model and the analytical model. While the analytical model is completely fixed at the end, the BCs of the FE model are a little more relaxed as they are imposed through cohesive element at only half the circumference (as opposed to a complete fixation of the entire cross-section). This may explain some of the differences seen between the two models.
The analytical model has a clear advantage over the FE model when it comes to computational cost. While the FE model required approximately one day of CPU time on a work station, the analytical model completed in a few seconds. Furthermore, the analytical model was capable of providing predictions for opening displacements beyond the point where the equilibrium iterations failed for the FE model. It may be possible to establish more efficient FE models than the one used in the present study, but it is our belief that the element length must be much smaller than the fiber diameter in order to suffi-ciently resolve the fracture process zone between the fiber and the matrix. Then the number of elements, and the computational cost, will be high even if beam elements are used. The authors have not been able to obtain good results with beam elements (where an underlying assumption is that the length of the element is much greater than the height and width), but this may be possible.
The finite element model predicts that the fibers in compression have a significant post-buckling capacity (see Fig. 9), while all forces are set equal to zero in the proposed analytical model when the critical load according to linearized buckling theory is reached. In this respect, the drop in tractions predicted by the proposed model is probably too abrupt, but on the conservative side.
As seen in Fig. 9(d), the finite element model predicts non-zero shear force in the initial phase for pure mode II. This is caused by the eccentricity of the axial loading of the fiber from the two rigid bodies some vertical distance apart (see Fig. 5). The effect of this is balanced out when there is an equal number of fibers bridging the two diagonal directions. Further, the finite element model doesn't predict   the axial force to instantly reach a plateau like the analytical model. This is because the fiber is modelled with finite thickness and stiffness, and is therefore allowed to shear between the two fracture surfaces. A more advanced beam model would not improve the predictions, as the bridging ligament is not acting as a beam at this stage, but rather as a bulky 3D-structure (the length is smaller than the height and width). Similarly, the finite element model predicts a horizontal force to be transferred to the fracture surface for very small mode I openings, see Fig. 9(a). This is due to shear stress at the interface between the fiber and matrix. Including shear deformations in the micromechanical model would not improve the predictions since the transferred force is assumed to be equal to the axial force acting at the neutral axis of the fiber. The error from this becomes negligible for larger openings. For an equal number of fibers bridging in the two directions, the effect would balance out also for small opening displacements.
Since bending and axial extension are the only admissible deformations of the micro-mechanical beam, the model predicts normal opening components to cause infinitely large normal tractions when the beam length approaches zero (see Fig. 9 and 10). This is a weakness of the proposed model. However, the model is not intended to predict the onset of fiber bridging at the crack tip but rather what happens in the wake of the crack tip where the fibers are of a finite length.
It should be noted that on the side of a circular fiber peeling off, the debond crack opening includes a mode III component and the opening mode will be a varying mixity of mode I and III along the circumference of the fiber. The results presented in the previous section were based on the assumption that the fracture energy is the same for all three modes, G Ic ¼ G IIc ¼ G IIIc ¼ G c . These are likely to be different in reality. Fracture energies dependent on mode mixity have not been implemented in the proposed model.
The model can be used to predict bridging forces from bundles of fibers being pulled out of the matrix instead of individual fibers. Appropriate values for EI and EA must then be used. However, adapting the fiber failure criterion is not straight forward.
The fibers are assumed to be prevented from rotations at the end of the free span. In reality, the two fracture surfaces may rotate relative to each other. The effect of this has been neglected. It should also be noted that the fibers will bridge between two points that did not initially coincide.
The effect of frictional sliding and interaction between fibers is neglected. Contact between crossing fibres may lead to point forces acting in their free span. This would affect the forces transferred between the two fracture surfaces through these fibres.
The low computational cost of the proposed model makes it well suited for parametric studies including stochastic variables. A small example of a parameter study has been shown in Section 4, but this is something the authors intend to study in more depth in the future. For pure mode II delamination, Eqs. (25)- (27) all give the same fraction of failed fibers. The analytical expressions for dissipated energy derived in Section 4 should therefore be exact. When there is a mode I component, on the other hand, Eq. (27) overestimates the fraction of failed fibers. The magnitude of the effects seen in Fig. 13 are therefore likely to be overestimated for mode I. Even if they may be quantitatively inaccurate, they qualitatively similar to those in pure mode II which should be correct.
An important finding in the parameter study was that increasing G c leads to more fibres failing earlier and a more rapid decrease in tractions. This caused reduced dissipated energy for both mode I and II in the short parameter study. To the best of our knowledge, this has not been shown analytically before.
The individual fibers are considered to only have stress/strain in the longitudinal direction (caused by bending, stretching or a combination thereof). No assumption regarding plane stress/strain in the matrix has been made since this is considered rigid in the model. Assuming a different width would therefore not affect the model predictions. In a real DCB test, for instance, the width would naturally affect the results, but that is at a different length scale than what is modelled here.
An important point is that δ x and δ y are only related to local opening displacements if the two fracture surfaces do not rotate relative to each other or deform within the span that the fiber bridges. Consider two rigid surfaces that are both prevented from rotating. Regardless of which points that are chosen on the two surfaces as reference points, the relative displacement will be the same for any arbitrary combination of tangential and normal movement.

Conclusions
A novel micromechanical model of crossover fiber bridging has been developed. The predicted macroscopic bridging laws of the new model were found to follow those predicted by a detailed FE model capable of taking large deflections and shear deformations into account.
It has been demonstrated that large deflections must be taken into account when modeling cross-over fiber bridging. Models based on linear beam theory will not be able to predict the correct stress and strain inside the fiber.
The proposed model allows strains to be integrated over the exposed surface of the fiber, and thus enables fiber fracture to be handled as a stochastic and size dependent process.

Data Availability
No data was used for the research described in the article.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.