Calcium phosphate cement reinforced with poly (vinyl alcohol) fibers An experimental and numerical failure analysis

Calcium phosphate cements (CPCs) have been widely used during the past decades as biocompatible bone substitution in maxillofacial, oral and orthopedic surgery. CPCs are injectable and are chemically resemblant to the mineral phase of native bone. Nevertheless, their low fracture toughness and high brittleness reduce their clinical applicability to weakly loaded bones. Reinforcement of CPC matrix with polymeric ﬁbers can overcome these mechanical drawbacks and signiﬁcantly enhance their toughness and strength. Such ﬁber-reinforced calcium phosphate cements (FRCPCs) have the potential to act as advanced bone substitute in load-bearing anatomical sites. This work achieves integrated experimental and numerical characterization of the mechanical properties of FRCPCs under bending and tensile loading. To this end, a 3-D numerical gradient enhanced damage model combined with a dimensionally-reduced ﬁber model are employed to develop a computational model for material characterization and to simulate the failure process of ﬁber-reinforced CPC matrix based on experimental data. In addition, an advanced interfacial constitutive law, derived from micromechanical pull-out tests, is used to represent the interaction between the polymeric ﬁber and CPC matrix. The presented computational model is successfully validated with the experimental results and offers a ﬁrm basis for further investigations on the development of numerical and experimental analysis of ﬁber-reinforced bone cements.

Reinforcement of brittle matrices with high-strength fibers has been extensively investigated for different applications in civil, aerospace and mechanical engineering. Due to their superior properties compared to conventional structural materials, fiberreinforced composites (FRCs) have attracted much attention in the past decades [36][37][38][39][40] . Previously, it was shown that the fiber reinforcement technique can also be effectively employed to enhance the mechanical properties of the cements used in orthopedic and dental applications [41][42][43][44][45][46][47][48][49][50] . However, fiber-reinforced calcium phosphate cements (FRCPCs) are still poorly investigated and understood in terms of the mechanism by which polymeric fibers reinforce the CPC matrix. Combined numerical and experimental studies describing the mechanical performance of FRCPCs under clinically relevant bending and direct tensile testing are not yet available. Hence, the main goal of this study is to develop an advanced, combined experimental-numerical model to investigate the mechanical behavior of fiber-reinforced CPCs under bending and direct tensile loading. To this end, the three main phases of FRCPCs should be considered, i.e. i) the fiber-matrix interface, ii) the CPC matrix, and iii) the dispersed fibers. Previously, we presented a combined experimental-numerical approach to describe the affinity between poly (vinyl alcohol) (PVA) fibers and the CPC matrix [51] . PVA fibers were selected in view of their highly effective reinforcing efficacy in civil engineering as well as the favorable biological performance of poly (vinyl alcohol) as biomaterial [52] . In this study, experimental micro-mechanical pull-out tests were combined with a numerical finite element (FE) model including distinct representation of the fiber, matrix and fiber-matrix interface with a predictive interfacial constitutive law [51] . The proposed interfacial constitutive law was validated experimentally and enabled prediction of all three main phases of the experimentally observed pull-out response, i.e. the elastic, debonding and frictional pull-out phases. Subsequently, we investigated the failure behavior of fiber-free calcium phosphate cements under bending and tensile loading by combining experimental tests and numerical modeling [53] . We employed a three dimensional (3-D) gradientenhanced damage model to computationally model the failure behavior of fiber-free CPC matrices in a mesh-objective and accurate manner [54] . This current study is the culmination of our two previously published experimental-numerical studies on the fibermatrix interface and the CPC matrix, respectively [51,53] . Based on these studies, we now present a fully integrated experimentalnumerical characterization of the mechanical performance of FR-CPCs. To this end, we use these two previously developed models combined with a dimensionally-reduced fiber model [55] to construct a complete computational framework for material characterization of FRCPCs.
In the presented model, the gradient-enhanced damage model captures the process of crack propagation in the fiber-free CPC matrix [53] . PVA fibers are assumed as linear elastic material, where their slippage is considered only in the axial direction and modeled using the dimensionally-reduced fiber model to allow a higher computational cost reduction [55,56] . As described in our previous work, a very fine mesh discretization is needed to represent the material length scale of the fiber-free CPC matrix and the corresponding simulations are numerically costly [53] . Therefore, we used the dimensionally-reduced fiber model in which fibers and matrix discretization are independent from one another [55] , inspired by the embedded reinforcement technique [57][58][59] , and the fibers are assumed as one-dimensional objects. In comparison to 3-D finite element models of fiber-reinforced composites where the exact geometry of fibers should be discretized, one-dimensional fiber models drastically reduce computational complexity and time of the numerical simulations.
Herein, we present the use of above-described models, i.e. i) the gradient-enhanced damage model, ii) three-phase fiber-matrix interface model and, iii) the dimensionally-reduced fiber model as an efficient computational tool to study the mechanical performance and failure behavior of FRCPCs. Moreover, we validate the presented computational model against experimental data obtained by subjecting FRCPCs to a range of three-point bending and direct tensile tests. Overall, this research offers a perspectives for the design and optimization of FRCPCs to mature their development as functional and reliable load-bearing biomaterials.

Three-point bending and tensile tests
Poly (vinyl alcohol) (PVA) fibers were utilized to reinforce CPC matrices in order to evaluate how fiber-matrix interfacial properties affect the macro-mechanical properties of fiber-reinforced CPCs. To perform the experimental three-point bending and direct tensile test, PVA fiber-reinforced CPC specimens were fabricated by first mixing 98 . 5 % α-tricalcium phosphate ( α-TCP, D50 of 2 . 97 μm , D90 of 6 . 06 μm , volume mean diameter of 3 . 51 μm ) (CAM Bioceramics, Leiden, the Netherlands) with 2 . 5 wt % PVA fibers (Kuraray Europe GmbH, RFS 400 / 18 mm , Hattersheim am Main, Germany, 200 μm in diameter and 18 mm in length) that were manually cut into lengths of 8 and 4 mm . Subsequently, a 4 wt % NaH 2 PO 4 . 2H 2 O (Merck, Darmstadt, Germany) aqueous solution was added at a liquid-to-powder ratio (L/P) of 1:2 and fully mixed for 1 min until a cementitious paste with randomly dispersed PVA fibers was formed [60,61] . The initial and final setting times of these cement formulation were ∼ 3 and ∼20 min, respectively [62,63] . The paste was subsequently cast into polydimethylsiloxane (PDMS) rectangular molds (40 × 10 × 10 mm ) , clamped between two glass slides and allowed to set at room temperature for 24 h . Afterwards, the specimens were removed from the molds and incubated in a phosphate-buffered saline solution (PBS) at 37 • C for 72 h to allow the CPC to fully cure. The weight of each coated PVA fiber ( l f = 18 mm ) was 0 . 77 mg , and each fiberreinforced CPC sample contained 125 mg of PVA fibers. The density of PVA fiber was ∼ 1 . 3 g / cm 3 and the fiber content was kept constant at 2 . 5 wt % for both the three-point and tensile tests irrespective of the fiber length. As a result, for both the three-point bending and tensile tests, the number of PVA fibers within the each fiber-reinforced CPC sample was calculated as 370 and 740 for fiber lengths 8 and 4 mm , respectively.
To ensure a controlled crack initiation and propagation process during the mechanical tests, a 3 × 1 mm single edge notch and two opposing 2 . 5 × 1 mm notches were cut into the specimens using a diamond-tipped circular saw blade. These notches prevent an uncontrolled and random failure process and allows the validation of the computational model parameters according to the experimental data. To ensure that all experimental tests were performed under hydrated conditions, the specimens were stored in 37 • C (PBS) until the moment of testing. Schematic sketches of both three-point bending and tensile test setups, boundary conditions and specimen geometry are illustrated in Figure SM Using a universal testing machine (LLOYD material testing, LS1 series) equipped with a 10 0 0 N load cell, the three-point bending test and the tensile test were performed at a crosshead speed of 1 mm / min . To perform the three-point bending tests the specimen was placed on two supporting pins with a span length of 30 mm and the load was applied incrementally at the mid-span. Values of fracture toughness, K IC , were computed using Eq. (1) and (2) based on ASTM C 1421 − 10 (an established test method for fracture toughness of advanced ceramics at ambient temperature [64] ) as: where F max is the maximum force and s 0 and b represent the support span and the specimen thickness, respectively. The notch depth is equal to a, w is the width of the test specimen and g is a function of ratio a/w, as shown in Eq. (2) . The bending strength f f l was calculated using Eq. (3) : where w denotes the distance between the tip of the notch and the top of the specimen. For tensile tests, the double-notched specimens was first mounted to the clamps of the universal testing machine. To this end, two plastic T-shaped bars were 3-D printed and glued to either end of the specimen using a two-component Pleximon glue (Evonik Röhm GmbH, Darmstadt, Germany). The specimen was subsequently placed into the testing machine and subjected to axial tensile loading until rupture. Here the tensile strength f ts was calculated as: where w is the distance between the tip of the two notches for the tensile test specimens. The work of fracture under bending ( WOF b ) and tensile ( WOF t ) loading were measured by dividing the total area under the force-displacement curves by the crosssectional surface area [61] .
The experimental data were reported as average ± standard deviation and analyzed statistically by means of one-way analysis of variance (ANOVA) followed by a Tukey post hoc test. For all tests, force-displacement curves were recorded for at least 10 samples per experimental group ( n ≥ 10 ) and a data collection frequency of 16 kHz was applied to record the experimental data points. The spatial distribution of fibers within test specimens was characterized using nano-computed tomography (nano-CT, Phoenix NanoTom M, General Electric, Wunstorf, Germany). Nano-CT analysis was obtained using a voxel size of 5 . 6 μm , vocal size spot of 0 . 84 mm , X-ray source of 70 kV / 200 μA , and exposure time of 500 ms without the application of a filter.

Results and discussion
For the three-point bending test and the tensile test the applied force was recorded as a function of the vertical displacement of the top-face of the specimen at mid-span section (bending tests) and the axial displacement of the top-end of the specimen (tensile tests). The corresponding force-displacement curves for threepoint bending and tensile tests of specimens reinforced with 8 and 4 mm PVA fibers are illustrated in Figs. 1 and 2 .
Generally, for both bending and tensile force-displacement curves of fiber-reinforced CPCs, three main phases can be distinguished. In the first phase, the specimen remains in its elastic regime. The incorporation of PVA fibers enhance the material resistance to deformation. In the second phase, by increasing the force, nano-and micro cracks form near the notch corners where the highest stress concentration occurs. Subsequently, the accumulation of these nano-and micro-cracks leads to a single macro-crack. In the third phase, toughening mechanisms become active and, in the crack wake, PVA fibers partially or fully bridge the crack to stabilize the crack propagation process. The bridging fibers debond and transmit the force across a crack and at the fiber-matrix interface. When further crack propagation occurs, fiber pull-out results in additional frictional sliding resistance, fibrillation of fiber surface, and finally slip-hardening behavior. Fiber type, number, distribution, embedded (debonded) length and location relative to the damage zone significantly affect the post peak response. A typical force-displacement curve obtained during the three-point bending test of a fiber-reinforced CPC specimen containing 8 mm PVA fibers is depicted in Figure SM Compared to the specimens reinforced with 8 mm PVA fibers, a smaller amount of energy was dissipated during the fracture of the fiber-reinforced CPC specimens with 4 mm PVA fibers (see Figs. 1 and 2 ). In our previous work, using the micromechanical pull-out test, the critical embedded length (l c ) of PVA fibers with 200 μm diameter embedded in a CPC matrix was 8 mm [51] . This implies that PVA fibers with 8 mm embedded length are the longest fiber that can be pulled out from their surrounding CPC matrix without rupture. This results in higher values of pull-out work, a longer interfacial frictional resistance phase and thus more energy dissipation. Consequently, the fiber length is a prominent parameter that can strongly affect the behavior of the post-peak regime and the efficiency of the interfacial properties in fiber-reinforced CPCs. Moreover, the presence of the fibers affected the damage width and propagation pattern. Images of two representative samples after execution of the bending and tensile tests are shown in Fig. 3 . Compared to the rapid failure induced by a narrow single macro-crack as observed in fiber-free CPC tests, we observed a slower crack growth process and a wider damage width [53] . A summary of the obtained experimental results is reported in Table 1 . In this table l e and d f are the fiber embedded length and fiber diameter, respectively. The fracture toughness, bending strength, tensile strength, bending and tensile modulus are denoted as K IC , f f l , f ts , E f l and E ts respectively. The values of WOF b and WOF t represent the work of fracture under bending and tensile loading, respectively.
Fracture toughness (see Eq. (1) ) is a quantitative parameter that expresses the material's resistance to crack propagation. Compared to the previously reported value of fracture toughness for fiberfree CPC matrix of similar chemical composition ( 0 . 17 MPa m 1 / 2 ), PVA fibers were able to increase the fracture toughness of fiberreinforced CPC matrices during the crack initiation and growth around 9% for 4 mm fibers and 30% for 8 mm fibers [53] . This low-fold increase might be related to sub-optimal fiber dimensions (low fiber aspect ratio for effective reinforcement) which were necessary for this study to be able to perform combined experimental and numerical analysis. The values measured for the fracture toughness of FRCPC matrices fit within the range reported for trabecular bone ( 0 . 1 − 0 . 8 MPa m 1 / 2 ), although they are still below values reported for cortical bone ( 2 − 12 MPa m 1 / 2 ). Parameters such as the interface affinity, fiber length, fiber volume fraction, fiber alignment and distribution directly affect the mechanical properties of FRCPCs [47,65,66] . Tailoring these parameters using computational models as developed herein will contribute to optimization of the mechanical properties of fiber-reinforced CPC matrices, which might ultimately lead to biomaterials with similar properties as found in bone.
Moreover, in comparison to fiber-free CPC specimens, the work of fracture values were significantly enhanced [47,53] . This means that more energy was dissipated during the fracture process of fiber-reinforced specimens, due to the crack-bridging and fiber pull-out process.
The nano-CT images of specimens reinforced with 8 and 4 mm PVA fibers are presented in Fig. 4 . The images show various sideviews of different cross sections of specimens after executing the mechanical tests. The light gray, dark gray and dark red colors represent the CPC matrix, PVA fibers and fracture surface, respectively. The nano-and micro-pores were homogeneously distributed over the specimens. PVA fibers were randomly distributed with respect to locations and orientations. As depicted in Fig. 4 -b and c, the PVA  fibers bridged the macro-cracks and stopped the crack propagation process. In addition, as shown in Fig. 4 -a and e, some of the PVA fibers near the damage zone were partially pulled out from the matrix followed by continuation of crack propagation. By comparing the fiber length of (partially) pulled out fibers in the nano-CT images and experimental specimens, it was concluded that none of the PVA fibers ruptured during the fracture process and the macrocrack mainly grew through the CPC matrix. This is in agreement with our micro-mechanical pull-out study where a critical embedded length of 8 mm was observed for PVA fibers [51] . Hence, the CPC matrix is considered as a homogeneous material, while the PVA fibers are assumed as linear elastic material in the numerical studies described in Section 3 .
To study the reliability of the measured material strength, Weibull modulus values were calculated for each experimental group [67,68] . The analysis was performed according to the ASTM standard C1239 − 07 using the following equation [69] : in which where P f and m are the failure probability at stress σ and the Weibull modulus, respectively. The value σ 0 is the characteristic stress at which 63% of the specimens fail [70] . This value for specimens containing 8 and 4 mm embedded fibers was measured as 4.6 and 3 . 6 MPa for bending test, and 1.5 and 1 . 2 MPa for the tensile test, respectively. The parameters n and i denote the total number of experimental specimens per experimental group and the specimen rank in ascending order of failure stresses. The Weibull plots of the bending and tensile strength of CPC matrix reinforced with 8 and 4 mm PVA fibers are presented in Fig. 5 . A summary of Weibull parameters are listed in Table 2 . The correlation coefficient is presented as R 2 .
As shown in Fig. 5 , the strength values for both three-point bending test and tensile test are approximately Weibull distributed.
The measured Weibull moduli (m ) for reinforced CPCs are higher compared to the Weibull moduli (m ) for fiber-free CPCs (9.7 for bending and 3.6 for tensile tests), which indicates that the reliability of fiber-reinforced samples was higher [53] . Moreover, the reinforced CPCs containing 4 mm PVA fiber displayed a higher Weibull modulus compared to CPCs reinforced containing 8 mm PVA fibers. This phenomenon can be explained by the fact that the larger number of smaller fibers were more uniformly distributed throughout the matrix and create a more homogeneous material.

Numerical modeling
In this section a gradient-enhanced damage model combined with a dimensionally-reduced fiber model is described to characterize the failure response of FRCPC under bending and tensile loading using finite element modeling. Herein, the main goal is to calibrate the numerical model parameters according to the obtained experimental results and present a numerical model that can be efficiently used to fully characterize the mechanical performance of FRCPCs.

Constitutive relations
We employ the 3-D gradient-enhanced damage model to numerically represent the properties of CPC matrix, while the modeling of PVA fibers was performed by means of a dimensionallyreduced fiber model. Furthermore, a three-phase traction separation law (TSL) is employed to connect the matrix and fiber elements and represents the fiber-matrix interface behaviour [51] .

CPC matrix
To numerically model the properties of the CPC matrix, we used the implicit gradient-enhanced damage model that was developed by Peerlings et al. [54] . This model solves the pathological mesh dependence, which is inherently present in standard local damage models. Herein, for small deformations, the rate form of stressstrain relation for the CPC matrix reads as: where σ m and ε m are stress and strain tensor and H m and D represent the fourth-order elasticity tensor and damage scalar variable, respectively. The damage parameter D varies between 0 for fully intact material and 1 for fully failed material. In addition to the equilibrium equation, the governing field equations include the diffusion equation (modified Helmholtz equation) [71] : where ∇ is the gradient operator and ē and ˜ e denote the nonlocal and local equivalent strain, respectively. The constant c is the material parameter (gradient parameter) of the dimension material length scale ( l m ) square. The natural boundary condition is used as follows: where n represents the unit outward normal. The damage parameter D in Eq. (7) is an explicit function of history parameter κ. The initial threshold of history parameter is κ i and its evaluation according to the Kuhn-Tucker relations is defined as [72,73] : The exponential damage evolution law is used to define the damage behavior [74,75] : in which α and β are the material parameters. The parameter β corresponds to the rate of damage growth. In order to determine the equivalent strain in Eq. (8) , the modified von-Mises definition was employed [54] : where I 1 and J 2 are the first and second invariant of the strain tensor and deviatoric strain tensor.

PVA fiber
According to the nano-CT observations, PVA fibers can be considered as linear elastic, while the effect of fiber surface fibrillation during the pull-out process was assumed negligible and was therefore ignored [51] . This can be also explained by the fact that the fiber length does not exceed the fiber critical embedded length ( l e ≤ l c ) and therefore the pull-out process occurs without fiber rupture. Hence, the stress-strain relation for the PVA fibers is isotropically linear elastic and can be generally defined as: where H f is the elastic stiffness tensor of a PVA fiber and σ f and ε f are the stress and strain tensor.

Interface between PVA fiber and CPC matrix
To introduce the interfacial properties between the PVA fibers and CPC matrix, we employed zero-thickness interface elements with a constitutive traction separation law (TSL). The interfacial traction t c along the fiber is a function of the displacement jump: where t c = [t s , t t , t n ] and T and [[ ˙ u int ]] are the cohesive tangent matrix and the rate form of displacement jump in local coordinate system attached to the interface surface. The cohesive tangent matrix T can be written in matrix form as: where T s , T t and T n are the shear and normal stiffnesses of the interface. Herein, we only consider fiber slip along the fiber axis direction, thereby the values of T t and T n , stiffnesses of the interface in the direction normal and perpendicular to the interface surface and fiber axis, are set artificially high to prevent interface separation in those directions. The displacement jump denote the slip between fiber and matrix s in tangential direction. The precise form of the interface tangential stiffness T s = ∂ t s /∂ s depends on the TSL used. In our previous work [51] , we performed micromechanical pull-out experiments to investigate the affinity between the PVA fibers and the CPC matrix. As a result, we proposed the use of a three-phase TSL to numerically model the complete pull-out process i.e. the elastic, debonding and frictional pull-out phases. According to the three-phase TSL formulation, the bond stress-slip relation reads as: (16) in which G represents the corresponding relative bond modulus and γ 0 , γ 1 , and γ 2 denote the parameters controlling the descending and ascending branches in the debonding and frictional stages. The starting point of the debonding phase and the sliding phase are presented as s 0 and s 1 , respectively (see [51] for more details).

Weak form of governing equations
The weak forms of equilibrium and diffusion equations governing the deformation and failure processes of a fiber-reinforced calcium phosphate cement specimen are presented below. The principle of virtual work in a domain with the boundaries of the domain can be described as follows (see Fig. 6 (21) in which, δu is the virtual displacement applied over the Dirichlet  [55,76] . The weak form of the diffusion equation (see Eq. (8) ) can be formulated as: Eqs. (17) and (22) form the governing system of coupled set of partial differential equations. In the following, the general form of the global system of equations for the finite element implementation is described.

Spatial discretization
Following the Galerkin approach, the displacement vector for matrix u m and fiber u f and the nonlocal equivalent strain ē can be discretized by means of shape functions as: 23) where matrix N m , N f and row vector N e contain the shape functions for the displacement fields and the non-local equivalent strains. Moreover, u ˜ m , u ˜ f and ē ˜ represent the nodal degrees of freedom of the matrix and fiber displacement components and the nonlocal equivalent strains, respectively. For each interface element the displacement jump can be determined as:  (26) where N f , 1 and N f , 2 are the shape functions of a one-dimensional fiber element (node 1 and 2) and N m , n is the shape function at node n of the parent matrix element [55] . Differentiation of Eq. (23) leads to the strain components and gradient of non-local equivalent strain: where B m , B f and B e contain the derivatives of the shape functions N m , N f and N e , respectively.

Linearization and the incremental-iterative solution procedure
The next step consists of the linearization of the governing equations in order to construct a consistent incremental-iterative Newton-Raphson solution procedure. In this regard, linearization at iteration i + 1 at the nodal level leads to: and at integration point level we have: in which, where represents the iterative increment between iteration i and i + 1 . In Eq. 30 , the relation between the history parameter κ and ē ( ∂ κ/∂ ē ) is equal to 1 when ē > κ 0 and 0 otherwise. Herein, κ 0 denotes the converged value of the history parameter in the previous increment. The set of governing Eqs. (17) and (22) at iteration i + 1 can be expressed as: where f ext and f int denote the discrete balance of external and internal nodal forces, respectively. Similar expressions are used for the f ext , e , f int , e in the diffusion equation for clarity ( f ext , e = 0 ). Substitution of expressions 28-30 into Eqs. 31 and (32) and following the standard linearization procedure gives (more details can be found in [54,55,77] ): The last step to construct the global system of equations, is introducing the stiffness matrices to rephrase the set of governing equations into a compact matrix form.

Matrix
Considering the terms related to the CPC matrix in Eqs. (33) and (34) , the following stiffness matrices can be introduced [72] :

Fiber
According to the dimensionally-reduced fiber model [55,56] , all PVA fibers were considered as one dimensional objects that can transfer the load only in axial direction. Fiber coupling is ignored for simplicity. In addition, since the fibers are discretized with 1-D elements, overlap of fibers is not considered. This idealization of fibers helps to drastically reduce the computational cost during the numerical simulations compared to the conformal finite element approach, however, still six extra degree of freedom per fiber element are added. The validity of the dimensionally-reduced fiber model compared to a conforming 3-D fiber mesh is discussed in detail in [55] . Using the dimensionally-reduced fiber model, the PVA fiber discretization is independent from the CPC matrix discretization and, moreover, there are no constraints on the number of PVA fibers that intersect a CPC matrix element. This simplifies the mesh discretization process. The PVA fibers are discretized using typical two-noded 3-D truss elements in which each node has three global degrees of freedom. To simplify the calculation process of the stiffness matrix for matrix elements that are crossed by fibers, the matrix contribution is calculated over the total volume of the element. Therefore an effective Young's modulus is used to generate the fiber stiffness matrix and cancel out the already computed matrix contribution in the fiber region. The stiffness matrix for each fiber element is given by [78] : (39) where E f , E m , A f and L e f denote the fiber and matrix Young's moduli, fiber cross section area and fiber element length, respectively. Using the node coordinates, the cosines can be defined as: (40) in which the global coordinate of fiber element ends are defined as ( x 1 , y 1 , z 1 ) and ( x 2 , y 2 , z 2 ).

Fiber-matrix interface
The cohesive tangent stiffness matrix for a given interface element between PVA fiber and CPC matrix (see Eqs. (24) and (33) ) can be formulated as [55,77] : and Considering fibers as one-dimensional objects, for each submatrix in Eq. (41) , the integration was computed over the length of interface element L e f (equal to fiber element length) and multiplied by C f that denotes the fiber cross-section circumference.

Global system of equations
By assembling all sub-matrices introduced in previous sections, the global structure of coupled system of equations ( Eqs. (33) and (34) ) can be written as: where the contribution of interface stiffness is shared between fiber and matrix degree of freedom. Eq. (47) is solved using an implicit Newton-Raphson scheme. In each iteration, the system of equations around the approximate solution calculated in the previous iteration is linearized and solved to compute a new approximation of the solution. This procedure is repeated until convergence of the solution is reached. Herein, the iteration superscript i + 1 and i related to the implicit Newton-Raphson solver are omitted for clarity. The proposed numerical model for solving this nonlinear system of governing equations is coded in the researchoriented C++ programming toolkit named Jive [79] . In the next section, the results of the numerical simulations are presented and the model is validated against the experimental results.

Results and discussion
A 3-D FE model was developed to numerically study the failure process of fiber-reinforced calcium phosphate cements under three-point bending test and tensile test. A 3-D model represents the fiber distribution in the CPC discretization more realistically despite the fact that it is computationally more costly compared to a two-dimensional (2-D) numerical framework.
In all numerical simulations the CPC matrix was discretized by tetrahedral continuum elements and truss elements were used to discretize the PVA fibers. The PVA fibers remain elastic and will not break during the pull-out process. Moreover, the zero-thickness interface elements with a three-phase TSL were incorporated into the model to represent the interfacial properties between the PVA fibers and CPC matrix. The problem boundary conditions were enforced as illustrated in Figure SM-1. For the three-point bending test, the specimen is subjected to a prescribed displacement at mid-span, the left supporting pin is fixed and the right pin is constrained in y-and z-directions. For the tensile test, the specimen is subjected to axial tensile loading while the bottom edge of the specimen is fixed.
To mimic the effects of PVA fibers being randomly distributed during the experimental tests, the PVA fibers with uniform fiber density distribution are arbitrary distributed in the CPC matrix discretization. For each numerical simulation, two random fiber distribution were considered. Furthermore, to achieve an accurate measurement of the stress profile over the vicinity of notches where the crack initiation and damage propagation are expected, we discretized this area with a finer mesh discretization compared to the rest of the structure.
The mesh discretization for both three-point bending and tensile tests of fiber-reinforced calcium phosphate cement specimens with two random PVA fiber distributions and two different fiber lengths ( l e = 8 and 4 mm ) and diameter d f = 200 μm are depicted in Figs. 7 and 8 , respectively, where h represents the element size in the refined region. In these figures, for clarity, just the outline and half of the matrix are illustrated for the first and second fibers distribution, respectively. Moreover, to show the differences, both fiber distributions, black and red lines, are depicted on top of each other (2-D view) for both bending and tensile simulations (see parts e and f in Figs. 7 and 8 ).
The numerical force-displacement curves for three-point bending and tensile tests of fiber-reinforced CPC with two random PVA fiber distributions and two different fiber embedded lengths ( l e = 8 and 4 mm ) together with the corresponding experimental results are plotted in Figs. 9 and 10 respectively. We selected the threepoint bending test for calibration purposes due to the more stable and controlled crack propagation during failure as compared to tensile testing [53] . In this regard the numerical model parameters were tuned in order to interpolate the average of maximum force and the displacement at peak load of three-point bending tests (black point in Figs. 9 -a and 10 -a) as closely as possible. The same set of calibrated parameters was used to predict the mechanical response of the fiber-reinforced CPC under the tensile tests. A B-spline interpolation technique was used to find a rep-   (12) ) [53] . The damage evolution variables α and β in Eq. (11) are set as 0.9 and 80, respectively. Moreover, the set of parameters for the interfacial traction separation law were employed according to the tuned parameters obtained from micro-mechanical pull-out tests presented in [51] .
The numerical results fit within the experimental envelope. This evidences that the proposed numerical model can reliably represent the bending and tensile responses of fiber-reinforced CPCs. The numerical results of both fiber distributions for each threepoint bending and tensile simulations are relatively similar. Small differences, mainly in the post peak regime, are related to the random fiber locations and therefore to a different level of fiber pull-    For both three-point bending and tensile tests of fiberreinforced CPCs, the damage profiles are wider compared to the fiber-free CPCs [53] . This was also observed experimentally by comparing the fracture surface of specimens after execution of the tests. For fiber-free CPCs, a narrow and localized single macrocrack with a smooth fracture surface was detected, whereas for fiber-reinforced CPCs a relatively wider and expanded crack with a tortuous fracture surface was observed (see also Fig. 3 ). In addition, the damage profile of the CPC matrix reinforced by PVA fibers with an embedded length of 8 mm is wider compared to those reinforced with a fiber embedded length of 4 mm (see Figs. 11 and 12 ). Due to the crack-bridging process, the crack can zip along the fibers and more fracture energy can be dissipated as a result of fiber pull-out. This can lead to a wider damage zone, larger values of work of fracture, and a more stable failure process.
Another parameter that can affect the width of the damage band is the material length scale, which is related to the average size of coalesced pores and material micro-structure imperfections and its effect was introduced by gradient parameter c in the diffusion equation to the system of governing equations. The effect of this parameter for CPC matrix is studied in our previous work [53] . In fact, the gradient parameter c in Eq. (8) adjusts the width of the damage band. Higher values of the gradient parameter correspond to a larger size of micro-structural voids and imperfections that lead to a wider localization zone. Furthermore, the material parameter β in the damage evolution law ( Eq. (11) ) determines the crack propagation rate and lower values of β results in a slower crack propagation and consequently a more ductile failure response. The improved mechanical properties of reinforced composites are attributed to the higher number of interfaces, thereby increasing the number of channels for crack propagation. In addition, the formation of additional microstructural defects produced by immature fiber-matrix interfaces facilitates crack deflection which leads to enhanced fracture energy dissipation and prevents the composite from catastrophic failure. Hence, to mimic the wider damage band and more stable damage growth that was experimentally observed for fiber-reinforced CPC, we selected a larger material length scale c and lower material parameter β for our simulations compared to our previous study on fiber-free CPCs [53] .
The good agreement between the numerical results and the experimental envelope evidences that the proposed combined numerical model can be reliably applied to further investigations on FRCPCs and replace the tedious experimental trial-and-error design procedures that are commonly used in biomaterials science and engineering.

Conclusions
In this work, the mechanical response of a CPC matrix reinforced by PVA fibers subjected to bending and tensile loading was investigated both experimentally and numerically. To develop a comprehensive numerical model, we proposed the use of 3-D gradient enhanced damage model combined with a dimensionally reduced fiber model to numerically characterize the failure behavior of a fiber-reinforced CPC matrix. The PVA fibers are idealized as one-dimensional objects which drastically reduced the computational cost during the numerical simulations compared to the conformal finite element approach. Moreover, an advanced traction constitutive law is employed to represent the fiber-matrix interfacial properties. We tuned the parameters of the combined numerical model according to the average of experimental data and showed that the proposed model is able to predict the failure response of a fiber-reinforced CPC matrix under bending and tensile loading with good accuracy. We showed that the PVA fiber embedded length is a key parameter that can affect the amount of fracture energy dissipation, damage zone and stability of failure process. We argue that the proposed combined numerical model can be employed with a good approximation for further studies on the development and optimization of fiber-reinforced CPCs.

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.