Experimental and numerical analysis on bending and tensile failure behavior of calcium phosphate cements

Since discovery injectable self-setting calcium phosphate cements (CPCs) are frequently used in orthopedic, oral and maxillofacial surgery due to their chemical resemblance to the mineral phase of native bone. However, these cements are very brittle, which complicates their application in load-bearing anatomical sites. Polymeric fibers can be used to transform brittle calcium phosphate cements into ductile and load-bearing biomaterials. To understand and optimize this process of fiber reinforcement, it is essential to characterize the mechanical properties of fiber-free calcium phosphate matrices in full detail. However, the mechanical performance of calcium phosphate cements is usually tested under compression only, whereas bending and tensile tests are hardly performed due to technical limitations. In addition, computational models describing failure behavior of calcium phosphate cements under these clinically more relevant loading scenarios have not yet been devel- oped. Here, we investigate the failure behavior of calcium phosphate cements under bending and tensile loading by combining, for the first time, experimental tests and numerical modeling. To this end, a 3-D gradient- enhanced damage model is developed in a finite element framework, and numerical results are correlated to experimental three-point bending and tensile tests to characterize the mechanical properties of calcium phos- phate cements in full detail. The presented computational model is successfully validated against experimental results and is able to predict the mechanical response of calcium phosphate cement under different types of loading with a unique set of parameters. This model offers a solid basis for further experimental and computational studies on the development of load-bearing bone cements.


Introduction
Calcium phosphate cements (CPCs) are highly valuable bone substitution materials that are frequently used in orthopedic, oral and maxillofacial surgery (Bohner et al., 2012;Zhang et al., 2014;Ginebra et al., 2010;Bohner, 2000;Habraken et al., 2016). Due to their excellent compatibility to bone tissue and their chemical similarity to the mineral content of bones and teeth, calcium phosphate cements have gained great interest in biomedical and clinical applications (Epple et al., 2010;Hench and Wilson, 1993;Sanzana et al., 2008). Consequently, CPCs have been used to treat e.g. craniofacial defects, bone fractures, periodontal defects and bone defects resulting from tumor resection (Friedman et al., 1998;Constantz et al., 1995;Al-Sanabani et al., 2013;Uchida et al., 1990).
CPCs are typically prepared at room temperature by mixing a precursor powder and an aqueous solution at a specific liquid-to-powder ratio (L/P ratio) (Thamaraiselvi and Rajeswari, 2004). Calcium phosphate cements set in situ at body temperature and form a paste that is moldable which can be directly injected into bone defects (Comuzzi et al., 2002). This procedure provides intimate adaptation to the irregular shape of bone defects. Consequently, CPCs represent attractive alternatives to conventional types of bioceramics which are difficult to shape and machine due to their non-self-setting nature . Although CPCs are widely used in various types of surgeries, these biomaterials are still associated with several drawbacks such as their poor mechanical performance, which currently hinders their clinical applications in load-bearing skeletal sites (Liu et al., 2013). The compressive strength of calcium phosphate cements typically ranges between 60 70 MPa, which exceeds values reported for native cancellous (trabecular) bone ð4 9 MPaÞ (Ikenaga et al., 1998;Ambard and Mueninghoff, 2006). Compressive, tension and torsion tests were performed to characterize the multiaxial mechanical properties of two type of CPCs, brushite and hydroxyapatite, by Charriere et al. (Charri� ere et al., 2001). They used a linear elasticity model and a Tsai-Wu failure criterion to describe the elastic and failure properties of CPCs. In their research the compressive strength of brushite and hydroxyapatite CPCs were reported as 10:7 and 75 MPa, respectively. Furthermore, the tensile strength of brushite and hydroxyapatite CPCs was measured as 1:3 and 3:5 MPa, respectively. Nonetheless, the applicability of CPCs in load-bearing sites is limited by their poor mechanical response under bending and tensile failure behavior (Dos Santos et al., 2000).
To improve the mechanical performance of CPCs, reinforcement using polymeric fibers is one of the most successful strategies. To this end, strong and stiff polymeric fibers are introduced into brittle CPC matrices to improve the fracture toughness of the material by arresting crack propagation (Krüger and Groll, 2012;Kucko et al., 2019). However, unlike fiber reinforcement of cementitious matrices in civil engineering, fundamental understanding of fiber reinforcement of calcium phosphate cements is still lacking. This also limits the efficient optimization of fiber reinforcement in terms of the properties, amount, dispersion and fiber-matrix affinity of the employed polymeric fibers.
The mechanical behavior of fiber-reinforced CPCs is determined by the three main phases, i.e. the CPC matrix, the dispersed fibers and the fiber-matrix interface. To understand and optimize this process of fiber reinforcement, it is essential to first obtain fundamental understanding of the mechanical properties of fiber-free calcium phosphate matrices under different types of loading . Diametral compression tests, also known as Brazilian tests, are frequently performed to measure the tensile strength of CPCs (Mirtchi et al., 1989;Bermudez et al., 1993;Ishikawa and Asaoka, 1995;Martin and Brown, 1995;Chow et al., 2000;Charriere et al., 2003;Luo et al., 2016). However, it should be realized this test is an indirect method to measure the tensile strength of brittle materials and therefore not the preferable testing method to determine these material characteristics (Krüger and Groll, 2012). Pittet and Lemaitre measured the tensile strength of CPC matrix using direct tensile and indirect diametral compression test (Pittet and Lemaitre, 2000). Using a Mohr circle approach, they showed that the diametral compression test generally produces values lower than the true tensile strength of the tested material. However, the mechanical performance of calcium phosphate cements is usually tested under diametral and uniaxial compression only, whereas clinically more relevant bending and direct tensile tests are hardly reported due to technical limitations.
Generally, tensile tests of highly brittle materials are complicated by several technical challenges which may reduce the accuracy of this test considerably. These technical challenges include axial misalignment, inappropriate load cell capacity, specimen imperfection, stress concentrations at fixtures, loading velocity and boundary conditions. These factors complicate the accuracy of tensile testing, especially for brittle materials with low tensile strengths . Furthermore, to ensure that specimens fail due to maximum tensile loading in the midsection, specific geometries are required such as dog-bone shaped or notched specimens, which are often difficult to fabricate using highly brittle materials. Bending tests are proposed as alternative, albeit with similar technical limitations for materials displaying highly brittle failure behavior. In addition, computational models describing failure behavior of calcium phosphate cements under clinically more relevant tensile or bending loading have not yet been developed.
Therefore, the main aim of the current work is to study the mechanical response of CPC matrices under tensile and bending loading using a combined experimental-numerical approach. To this end, we use the gradient-enhanced damage model to computationally model the failure behavior of CPC matrices (Peerlings et al., 1996a). The gradient-enhanced damage model has been used in civil engineering to predict material failure resulting from initiation and propagation of damage (De Borst et al., 1993Geers et al., 1999;Kuhl et al., 2000;Pham et al., 2011;Waffenschmidt et al., 2014;Poh and Swaddiwudhipong, 2009;de S� a et al., 2006;Pearce et al., 2004). Similar to concrete in civil engineering, failure in CPC matrices is caused by nucleation, propagation and coalescence of micro-cracks and microscopic defects that ultimately lead to dominant macroscopic cracks.
The numerical implementation of the gradient-enhanced damage model in the finite element method (FEM) is straightforward, and this so-called nonlocal approach is objective with respect to FE mesh discretization (De Borst et al., 1993). This means that with refining the finite element discretization the model converges to a solution with a finite energy dissipation and a finite damage bandwidth (Peerlings et al., 1996a). Hence, the use of the gradient-enhanced damage model is suggested as a genuine and robust numerical model to investigate failure behavior of calcium phosphate cements. The numerical model developed herein can provide a solid basis for the development of advanced numerical models to investigate and optimize the mechanical properties and failure behavior of fiber-reinforced CPCs.
In this work, we (i) develop an experimental protocol to perform bending and tensile tests for CPCs; (ii) calculate the failure probability of bending and tensile tests for CPCs, (iii) present the gradient-enhanced damage model as a robust and objective numerical model to study the mechanical performance and failure behavior of CPCs, and (iv) derive a unique set of numerical parameters to validate the numerical model against experimental data.

Three-point bending and tensile tests
To study the mechanical response of the CPC matrix under bending and tensile loading, test specimens were prepared using a previously established method by adding a 4 wt% NaH 2 PO 4 :2H 2 O (Merck, Darmstadt, Germany) aqueous solution to the CPC powder consisting of 100% α-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) . Subsequently, the two phases were thoroughly mixed at a liquid-to-powder ratio (L/P) of 1:2 for 1 min to form a cementitious paste. The initial and final setting times of cementitious paste were e3 and e20 min, respectively Castro et al., 2017;Zuo et al., 2010;Petre et al., 2019). Next, the paste was placed in a polydimethylsiloxane (PDMS) rectangular mold ð40 � 10 � 10 mmÞ, clamped between two glass slides and set at room temperature for 24 h. The rectangular test specimens were subsequently removed from the molds and incubated in a phosphate-buffered saline solution (PBS) at 37 ∘ C for 72 h to allow for full hardening of the CPC. In previous studies performed by our research group, the complete conversion of our CPC hydroxyapatite was confirmed using X-ray diffraction, as evidenced by the fact that all diffraction peaks corresponded to hydroxyapatite, whereas peaks characteristic of α-tricalcium phosphate were not detected (Kucko et al., 2019). In addition, the porosity was reported to range between 45 55 vol% based on both theoretical calculation and experimental characterization (Kucko et al., 2019;Lodoso-Torrecilla et al., 2018). Furthermore, the pore size distribution of this CPC was characterized using X-ray microtomography, which revealed that the main pore size was about 60 μm (Lodoso--Torrecilla et al., 2018).
To prepare specimens for three-point bending tests, a 3 � 1 mm single edge notch was cut into the specimens using a diamond-tipped circular saw blade. For tensile tests, two opposing 2:5 � 1 mm notches were cut on both sides of the specimens. These notches were applied since the failure process for unnotched specimens is unpredictable, resulting into failure of the CPC matrix at random locations induced by e.g. material inhomogeneity and/or stress concentrations at fixtures of the test setup. The presence of the notches guarantees a controlled crack propagation pattern, thereby allowing for validation of the numerical model parameters by means of experimental tests. The test specimens were analyzed using nanocomputed tomography (nano-CT, Phoenix NanoTom M, General Electric, Wunstorf, Germany). A nano-CT analysis was obtained using a voxel size of 5:6 μm, X-ray source of 70 kV= 200 μA, and exposure time of 500 ms without the application of a filter.
Representative sketches of both three-point bending and tensile test setups are depicted in Fig. 1. In order to ensure that all mechanical tests were performed under hydrated conditions, the specimens were stored in 37 ∘ C (PBS) until tests were performed.
Both the three-point bending test and the tensile test were performed using a universal testing machine (LLOYD material testing, LS1 series) equipped with a 1000 N load cell at a crosshead speed of 1 mm= min. For the three-point bending tests, the load was applied at mid-span of the specimen positioned with a support span of 30 mm. In order to perform tensile tests, the double-notched CPC specimens were first mounted to the clamps of the universal testing machine. However, mounting of the specimens directly to the clamps of the test setup was not feasible since this procedure generated excessive stress concentrations resulting into premature failure of the specimen. Hence, two plastic T-shaped bars were 3-D printed and glued (Pleximon) (Evonik R€ ohm GmbH, Darmstadt, Germany) to both ends of the specimens (see Fig. 1). Subsequently, the load was applied in the axial direction of the glued test specimens. For all tests a data collection frequency of 16 kHz was employed to record the experimental data points. Force-displacement curves were recorded for ten samples per experimental group (n ¼ 10).
Fracture toughness values (K IC ) were calculated according to ASTM C1421-10 (an established test method for fracture toughness of advanced ceramics at ambient temperature (ASTM, 1999)) as follows: here F max , s 0 and b denote the maximum force, support span and the specimen thickness. The notch depth is equal to a and w is the width of the test specimen.
The flexural strength f fl and tensile strength f ts were indirectly measured using Equations (3) and (4): where w' is the distance between the tip of the notch and the top of the specimen for bending test and w 00 is the distance between the tip of the two notches for the tensile test.
The experimental data were statistically analyzed by means of a oneway analysis of variance (ANOVA) followed by a Tukey post hoc test and are represented as average � standard deviation. All experimental data were analyzed for outliers using Tukey's criteria based on the maximum peak load.

Results and discussion
The experimental results for the three-point bending test and the tensile test are illustrated in Fig. 2. The applied force for bending tests was recorded against the vertical displacement (relative to top-face of test specimens at mid-span section). For the tensile test, the applied force was measured against the axial displacement of the points placed at the top-end of the specimen. The experimental machine was not able to trace the potential snap-back instability. As depicted in Fig. 2, both bending and tensile tests revealed a brittle failure response. The failure process started with initiation of micro-cracks around the corners of the notches where the highest stress concentrations occurred. Rapid propagation of these micro-cracks resulted in a single macro-crack leading to complete failure. Images of two representative samples after execution of the bending and tensile tests are shown in Fig. 3. Fig. 2 shows the typical scatter of experimental results for the bending test. The large variability of the obtained data is typical for testing of brittle ceramics and may be caused by several factors such as presence of material defects and inhomogeneities or misalignment of mechanical loading. However, the data variation was much larger for tensile tests than bending tests. Although it should be realized that tensile testing allows for reliable evaluation of tensile properties of ductile materials, tensile testing of brittle (bio)materials is extremely challenging. Tensile tests of brittle materials such as CPCs are highly sensitive to various experimental factors including misalignment of mechanical loading. For instance, mounting of CPC specimens to the tensile bench (see Fig. 1) could only be achieved practically by gluing printed plastic T-shaped bars to the specimens. However, these bars might have contributed to misaligned loading conditions, thereby increasing data variation. Therefore, we attribute the large data variation of the tensile tests to the effects of these T-shaped bars and Pleximon glue on load transfer to the CPC matrix.
A summary of experimental results is listed in Table 1. In this table F max , x max p and K IC represent the maximum peak load, maximum extension at peak load, and fracture toughness, respectively, while, the f fl and f ts correspond to the flexural and tensile strength, respectively. The values of E fl and E ts denote the flexural and tensile modulus, respectively.
The average value of fracture toughness is around 0:17 MPa m 1=2 , which corresponds to the range of 0:01 0:32 MPa m 1=2 as reported previously for CPC (Zhang et al., 2006(Zhang et al., , 2011Weir and Xu, 2010;Morgan et al., 1997;Harmata et al., 2015). Due to this very low fracture toughness, the reliability of CPC is compromised by the inherent presence of nano-and micro-sized pores. Although nano-and micropores can stop crack propagation and may lead to enhanced toughness of ductile materials, pores in brittle materials mainly act as crack nucleation sites causing uncontrolled and sudden failure. The average value of measured flexural modulus is in the range of previously reported values for different CPC formulations (0:4 8 GPa) Zhang et al., 2012;Boroujeni et al., 2013;Liu et al., 2014;Rajzer et al., 2016;. It should be noted that the mechanical properties of brittle materials as measured using pure tensile test are usually lower than the corresponding properties measured using bending test. During the tensile test the entire volume and cross section of the specimen is subjected to the maximum tensile stresses, while during the bending test the maximum tensile stresses are concentrated in a smaller volume under the neutral axis, which results in a lower probability of failure. The reported fracture toughness values for trabecular and cortical bone are around 0:1 0:8 MPa m 1=2 and 2 12 MPa m 1=2 , respectively (Fu et al., 2011). Consequently, the low fracture toughness of CPC relative to the bone tissue hampers the application of this biomaterial in load-bearing anatomical sites, and emphasizes the need for effective reinforcement strategies to develop advanced CPCs with enhanced toughness and strength (Dos Santos et al., 2000).
To better understand the statistical variation of the measured flexural and tensile strength, the Weibull modulus values were measured (Weibullet al., 1951;Fisher and Tippett, 1928). Weibull statistics are commonly employed to predict the failure probability distribution of brittle materials (Barralet et al., 2002a;Fleming et al., 2003;Dowling and Fleming, 2008). The Weibull modulus of porous calcium phosphate bioceramic scaffolds tested under compression was reported to range between 3 to 10 (Martínez-V� azquez et al., 2010; Bose et al., 2013), whereas a Weibull modulus of 5.6 was reported for apatite cements tested using four-point bending tests (Barralet et al., 2002b). Generally, however, the analysis of the Weibull modulus for CPCs has received little attention as compared to various other types of brittle materials. Herein, the Weibull modulus evaluation for three-point bending tests and tensile tests was performed according to ASTM standard C1239-07 using the following equation: where Fig. 2. Force-displacement curves for the three-point bending test (a) and the tensile test(b) tests of rectangular CPC specimens ð40 � 10 � 10 mmÞ.
Herein P f and m denote the failure probability at stress σ and the Weibull modulus, respectively. The value σ 0 is the characteristic stress at which the probability of failure is around 63% (Wachtman et al., 2009). The value n is equal to the total number of experimental specimens and i represents the specimen rank in ascending order of failure stresses. The Weibull plots of the measured flexural and tensile strength are illustrated in Fig. 4. A summary of Weibull parameters are listed in Table 2. The values of Weibull modulus (slope of the statistical trend-line in Fig. 4) for flexural and tensile strength were measured as 9.7 and 3.6, respectively.
The Weibull modulus of flexural strength was approximately 2.7 times larger than the Weibull modulus of tensile strength. Larger Weibull modulus values correspond to a more reliable and predictable failure process. This Weibull analysis confirms that the mechanical performance of our CPC specimens was much more reliable upon flexural testing as compared to tensile testing.
As shown in Fig. 5, pores of test specimens were more or less homogeneously distributed over the CPC matrix, and the size of the nanoand micro-pores was small compared to the macroscopic length scale of the specimens. Therefore, the CPC matrix was assumed as a homogeneous material in the numerical studies in Section 3.

Numerical modeling
In this section we propose the use of a gradient-enhanced damage model to study the failure behavior of CPC matrix under bending and tensile loading. The main goal is to uniquely characterize the numerical model parameters based on the obtained experimental data.

Gradient-enhanced damage model
The implicit gradient-enhanced damage model was developed by Peerlings et al. (1996a), where a discretized weak form of a system of coupled differential equations was solved using an incremental-iterative Newton-Raphson approach . Herein, the first equation in the system of coupled differential equations is the classical equilibrium relation that can be written in component form as: where σ ij and ε kl are the Cauchy stress and strain tensors, and C ijkl denotes the fourth-order elasticity tensor given by: in which λ and μ are the Lame's constants and δ denotes the Kronecker delta: where E and ν are the Young's modulus and Poisson's ratio of the material, respectively. The parameter D in Equation (7) is the isotropic damage variable that represents the effects of micro-structural defects in the material. The damage parameter can vary between 0 and 1, where D ¼ 0 corresponds to undamaged material and D ¼ 1 defines fully damaged material where the material has lost its stress-bearing capacity. Values of D between 0 and 1 correspond to crack propagation and development of material damage. The damage variable is an explicit function of history parameter D ¼ DðκÞ (with initial value κ ¼ κ i ). In the gradient-enhanced damage model, as a non-local damage model, the history parameter κ is related to the non-local equivalent strain ε eq and its evaluation can be mathematically defined by the Kuhn-Tucker relations Peerlings et al., 1998):

Fig. 5. Nano-CT of the cross section of a CPC matrix.
A. Paknahad et al. In non-local damage mechanics the damage growth in a material point depends not only on the local deformation in that point itself, but also on the average deformation in its direct vicinity. This non-local effect can be incorporated into a system of coupled differential equations (see also Equation (7)) as a diffusion equation: ε eq cr 2 ε eq ¼ ε eq (12) where r 2 is the Laplacian operator and ε eq and c are the local equivalent strain and the gradient parameter, respectively. The gradient parameter c is a material parameter that represents the average size of material micro-structural defects (micro-pores, micro-cracks) or simply the material length scale l m (Peerlings et al., 1996b).
In order to characterize the failure behavior of a material, different definitions of the damage evolution law can be employed in Equation (7) (Mazars and Pijaudier-Cabot, 1989;Geers et al., 1997). In this paper, we used the exponential damage evolution law as follows (Mazars and Pijaudier-Cabot, 1989;de Borst et al., 1998): where κ i is the damage initiation threshold and α and β are the material parameters that should be tuned based on experimental data. The exponential damage law facilitates the numerical model to remain stable until the end of the failure process where ε→∞ and the stress (see Equation (7)) goes to ð1 αÞC ijkl : ε kl . The parameter β defines the damage growth rate. High values of β lead to a more brittle material response with faster crack propagation. The modified von-Mises definition for equivalent strain ε eq (see Equation (12)) was used for all simulations (De Vree et al., 1995;Peerlings et al., 1996a): where I 1 is the first invariant of the strain tensor and J 2 is the second invariant of deviatoric strain tensor, given by: and k represents the ratio between the compressive and tensile strength.
The natural boundary condition is used as follows: r ! ε eq : n in which n ! is the unit outward normal (Peerlings et al., 1996a).
The nonlinear system of equations is solved by the implicit Newton-Raphson scheme and, the numerical model is coded in the researchoriented Cþþ programming toolkit named Jive (Lingen and Stroeven, 2002).

Results and discussion
A 3-D FE model was constructed to analyze the above-described three-point bending test and tensile test. Although 3-D modeling is relatively costly from a computational point of view compared to a 2-D approach, the presented 3-D model allows for extension to 3-D models for fiber-reinforced CPC with a realistic representation of the fiber distribution in the CPC discretization. The CPC matrix for all simulations was discretized by tetrahedral continuum elements and the boundary conditions were applied as shown in Fig. 1.
In order to obtain an accurate measurement of the stress profile over the area where the crack initiation and damage propagation are expected, we discretized the region in the direct vicinity of notches with a finer mesh compared to the rest of the structure. To perform the mesh refinement study and assure that the model is mesh-independent, a mesh discretization with three different element sizes was used in all simulations. The mesh discretization of the central part of the specimens for both three-point bending and tensile tests is depicted for three different element sizes in Figs. 6 and 7 where h represents the element size in the refined region.
The three-point bending test has been selected for calibration purposes due to the higher Weibull modulus and lower standard deviation compared to the tensile test data. The numerical parameters are calibrated with respect to the data obtained using experimental three-point bending to predict the tensile behavior of the CPC matrix with the same optimized and unique set of parameters. The numerical forcedisplacement curves are plotted in Fig. 8 together with the corresponding experimental results.
The model parameters are tuned in order to interpolate the average of maximum force and the displacement at peak load of three-point bending tests (red point in Fig. 8-a) as closely as possible. The damage threshold is set as κ i ¼ 0:0035 and the damage evolution variables α and β in Equation (13) are taken as 1 and 100, respectively. The numerical simulations are performed with gradient parameter c ¼ 0:01 and ν ¼ 0:2 and k ¼ 18 in the modified von-Mises definition (see Equation (14)).
We observed a good agreement between the numerical data and the experimental envelope, which evidences that the proposed numerical model is able to capture the bending and tensile behavior of the CPC matrix with good approximation.
The mesh refinement study for three different mesh sizes h ¼ 0:2; 0:1 and 0:05 mm (see Figs. 6 and 7) are depicted in Fig. 9. The numerical results for the three-point bending test and the tensile test converged for the smaller mesh sizes h ¼ 0:1 and 0:05 mm. This indicates that the numerical model is indeed mesh independent. The evolution of the damage profile for the three-point bending test and the tensile test are illustrated in Figs. 10 and 11.
The material length scale is related to the average size of coalesced pores and material micro-structure imperfections that can affect the width of macro-cracks. Figs. 12 and 13 show how different values of the gradient parameter c ¼ 0:01; 0:02 and 0:03 (related to material length scale) affect the damage profile for the three-point bending test and tensile test, respectively. Higher values of the gradient parameter result in a wider damage zone.
A lower value of c leads to a more brittle response characterized by a more narrow and localized damage zone that is a more realistic representation of the single macro-crack as observed experimentally (see Fig. 3). In fact, the gradient parameter c, the internal length scale of the gradient-enhanced model, sets the width of the localization zone. Lower values of gradient parameter c imply a smaller damage zone and the use of a very fine mesh specifically around the notches where damage takes place. This made the simulations numerically more time-consuming due to the large number of degrees of freedom (DOF) in the system of equations.
The good agreement between the numerical predictions and the experimental data evidences that the proposed numerical model can be reliably used to study the performance of CPCs without the tedious trial and error approaches that were commonly employed.

Conclusions
In this work, we experimentally studied the response of a CPC matrix under bending and tensile loading. We proposed an experimental protocol to perform the three-point bending and tensile test for a CPC matrix. We determined the Weibull modulus of the flexural and tensile strength for the CPC matrix. To derive a numerical model, we proposed the use of a gradient-enhanced damage model to numerically characterize the failure behavior of a CPC matrix. It was shown that the proposed finite element model is objective with respect to mesh discretization. We calibrated the numerical model parameters and the material length scale according to the experimental data and showed    Fig. 9. Mesh refinement study for three-point bending (a) and tensile (b) tests of rectangular CPC specimens ð40 � 10 � 10 mmÞ. that the model is able to predict the bending and tensile responses of the CPC matrix within the experimental envelope, using a unique set of parameters. To obtain a realistic representation of the narrow and localized single macro-crack, as observed experimentally, the value of the gradient parameter should be set sufficiently small which required a very fine mesh discretization and inevitably a large number of degrees of freedom. We argue that the proposed numerical model can be used with a good accuracy for further investigations on the development of CPCs.