Calcium phosphate cement reinforced with poly (vinyl alcohol) fibers

Delft University of Technology

Calcium phosphate cement reinforced with poly (vinyl alcohol) fibers

An experimental and numerical failure analysis

Paknahad, Ali; Goudarzi, Mohsen; Kucko, Nathan W.; Leeuwenburgh, Sander C.G.; Sluys, Lambertus J.

DOI

10.1016/j.actbio.2020.10.014

Publication date

2021

Document Version

Final published version

Published in

Acta Biomaterialia

Citation (APA)

Paknahad, A., Goudarzi, M., Kucko, N. W., Leeuwenburgh, S. C. G., & Sluys, L. J. (2021). Calcium

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

Acta Biomaterialia, 119, 458-471. https://doi.org/10.1016/j.actbio.2020.10.014

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

ContentslistsavailableatScienceDirect

Acta

Biomaterialia

journalhomepage:www.elsevier.com/locate/actbio

Calcium

phosphate

cement

reinforced

with

poly

(vinyl

alcohol)

fibers:

An

experimental

and

numerical

failure

analysis

Ali

Paknahad

a,b,∗

_,

_Mohsen

_Goudarzi

b

_,

_Nathan

_W.

_Kucko

a

_,

_Sander

_C.G.

_Leeuwenburgh

a

_,

Lambertus

J.

Sluys

b

a Department of Dentistry-Regenerative Biomaterials, Radboud Institute for Molecular Life Sciences, Radboud University Medical Center, Nijmegen, the Netherlands

b Faculty of Civil Engineering and Geosciences Delft University of Technology, Delft, the Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 3 July 2020 Revised 7 October 2020 Accepted 9 October 2020 Available online 24 October 2020 Keywords:

Fiber-reinforced calcium phosphate cements Three-point bending test

Tensile test Numerical modeling

a

b

s

t

r

a

c

t

Calciumphosphatecements(CPCs)havebeenwidelyusedduringthepastdecadesasbiocompatiblebone substitutioninmaxillofacial,oralandorthopedicsurgery.CPCsareinjectableandarechemically resem-blanttothemineralphaseofnativebone.Nevertheless,theirlowfracturetoughnessandhighbrittleness reducetheirclinicalapplicabilitytoweaklyloadedbones.ReinforcementofCPCmatrixwithpolymeric fiberscanovercomethesemechanicaldrawbacksandsignificantlyenhancetheirtoughnessandstrength. Suchfiber-reinforcedcalciumphosphatecements(FRCPCs)havethe potentialtoactas advancedbone substitute inload-bearinganatomicalsites. Thiswork achievesintegrated experimentaland numerical characterizationofthemechanicalpropertiesofFRCPCsunderbendingandtensileloading.Tothisend, a3-Dnumericalgradientenhanceddamagemodelcombinedwithadimensionally-reducedfibermodel areemployedtodevelopacomputationalmodelformaterialcharacterizationandtosimulatethefailure processoffiber-reinforcedCPCmatrixbasedonexperimentaldata.Inaddition,anadvancedinterfacial constitutivelaw,derivedfrom micromechanicalpull-out tests,is usedtorepresent theinteraction be-tweenthepolymericfiberandCPCmatrix.Thepresentedcomputationalmodelissuccessfullyvalidated withtheexperimentalresults andoffers afirmbasisforfurtherinvestigationsonthedevelopment of numericalandexperimentalanalysisoffiber-reinforcedbonecements.

© 2020ActaMaterialiaInc.PublishedbyElsevierLtd. ThisisanopenaccessarticleundertheCCBYlicense(http://creativecommons.org/licenses/by/4.0/)

Statement of significance

Reinforcement of calcium phosphate cements (CPCs) with polymericfibers,suchaspoly(vinylalcohol)(PVA)fibers,can substantially enhance their toughness and strength,thereby potentiallyovercomingthemechanicaldrawbacksofCPCs.In this study, we develop a combined experimental-numerical frameworktoinvestigatethemechanicalbehaviorandfailure process of fiber-reinforced calcium phosphate cements (FR-CPCs)underbendingandtensiletesting. Herein,wepresent afullyintegratedexperimental-numericalcharacterizationof the mechanical performance of FRCPCs considering three mainphases,i.e.i)thefiber-matrixinterface,ii)theCPC ma-trix,andiii)thedispersedfibers.Thepresentedmodelis

val-∗_{Corresponding author.}

E-mail addresses: a.paknahad@tudelft.nl , ali.paknahad@radboudumc.nl (A. Pak- nahad).

idatedexperimentallyandcanbe usedwithagood approxi-mationforfurtherstudiesonthedevelopmentand optimiza-tionofFRCPCs.

1. Introduction

Calciumphosphatecements(CPCs) arewidelyappliedin max-illofacial, oral and orthopedic surgery [1–6]. These cements are osteoconductive, biocompatible andchemically resemblant to the mineralphase of bonesandteeth [7–12]. Over thepast decades, injectable self-setting calcium phosphate cements have been in-creasinglyused to facilitate repair of various typesof bone frac-tures,periodontalandcraniofacialdefects [13–16].

However, despite their favorable biological performance, sev-eral crucial drawbacks are still associated to CPCs which limit their clinical applicability. In particular, the high brittleness and low fracture toughness of CPCs still limits their potential usage

https://doi.org/10.1016/j.actbio.2020.10.014

in highly loaded skeletal sites [17]. The compressive strength of CPCs is typically higher than native cancellous (trabecular) bone

(

1− 30MPa

)

andlower than cortical bone

(

95− 230 MPa

)

[18– 21]. Reported fracture toughness values of CPCs range between 0.01− 0.32 MPam1/2 _{[22–26]. Thesevalues are noticeably lower} than previouslyreportedfracturetoughness valuesofhuman cor-ticalbonethattypicallyrangebetween2_{− 12}MPam1/2 _[27]. Fur-thermore,CPCshavepoorresistancetotensileforces [2]. Nonethe-less, mechanicalpropertiesof CPCs wereassessed mostlyby per-forming diametral and uniaxial compression testsand their me-chanical response to bending and tensile loading are rarely re-portedduetotheinherenttechnicalchallenges [28–35].

Reinforcementofbrittlematriceswithhigh-strengthfibershas been extensively investigated for different applications in civil, aerospace and mechanical engineering. Due to their superior properties compared to conventional structural materials, fiber-reinforcedcomposites(FRCs)haveattractedmuchattentioninthe past decades [36–40]. Previously,it wasshown that thefiber re-inforcement technique can also be effectively employed to en-hance the mechanical properties of the cements used in ortho-pedic and dental applications [41–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 experi-mental studies describing themechanical performance of FRCPCs underclinicallyrelevantbendinganddirecttensiletestingarenot yet available.Hence, themain goalofthisstudyistodevelop an advanced, combinedexperimental-numerical modelto investigate the mechanical behavior of fiber-reinforced CPCs under bending and directtensile loading. Tothisend, the threemain phases of FRCPCs should be considered, i.e. i) the fiber-matrix interface, ii) the CPC matrix, andiii) the dispersed fibers. Previously, we pre-sented a combined experimental-numerical approach to describe theaffinity betweenpoly(vinyl alcohol)(PVA)fibers andtheCPC matrix [51].PVAfiberswereselectedinviewoftheirhighly effec-tivereinforcingefficacyincivilengineeringaswellasthefavorable biological performance ofpoly(vinyl alcohol)asbiomaterial [52]. In this study, experimental micro-mechanical pull-outtests were combined with a numerical finite element (FE) model including distinct representationof thefiber,matrix andfiber-matrix inter-face with a predictive interfacial constitutive law [51]. The pro-posedinterfacialconstitutivelawwasvalidatedexperimentallyand enabledpredictionofallthreemain phasesoftheexperimentally observed pull-out response, i.e. the elastic, debonding and fric-tional pull-out phases. Subsequently, we investigated the failure behavior offiber-free calcium phosphatecements underbending andtensileloadingbycombiningexperimentaltestsandnumerical modeling [53].Weemployed a threedimensional (3-D) gradient-enhanced damagemodeltocomputationallymodelthefailure be-havior of fiber-free CPC matrices in a mesh-objective and accu-ratemanner [54].Thiscurrentstudyistheculminationofourtwo previously publishedexperimental-numericalstudiesonthe fiber-matrix interface and the CPC matrix, respectively [51,53]. Based onthesestudies,wenowpresentafullyintegrated experimental-numerical characterization of the mechanicalperformance of FR-CPCs. Tothisend, weusethesetwo previouslydevelopedmodels combined witha dimensionally-reducedfiber model [55] to con-structacompletecomputationalframeworkformaterial character-izationofFRCPCs.

In thepresentedmodel,thegradient-enhanceddamagemodel capturestheprocessofcrackpropagationinthefiber-freeCPC ma-trix [53].PVAfibersareassumedaslinearelastic material,where their slippageisconsidered only in theaxial directionand mod-eledusingthedimensionally-reducedfibermodeltoallowahigher computational cost reduction [55,56]. As described in our previ-ous work, a very fine meshdiscretization is neededto represent

themateriallengthscaleofthefiber-freeCPCmatrixandthe cor-respondingsimulations are numericallycostly [53].Therefore,we used the dimensionally-reduced fibermodel in which fibers and matrix discretizationare independentfrom one another [55], in-spiredbytheembeddedreinforcementtechnique [57–59],andthe fibersareassumedasone-dimensionalobjects.Incomparisonto 3-Dfiniteelementmodels offiber-reinforcedcomposites wherethe exact geometry of fibers should be discretized, one-dimensional fibermodelsdrasticallyreducecomputationalcomplexityandtime ofthenumericalsimulations.

Herein,wepresenttheuseofabove-describedmodels,i.e.i)the gradient-enhanceddamagemodel,ii)three-phasefiber-matrix in-terfacemodeland,iii)thedimensionally-reducedfibermodelasan efficientcomputational toolto studythemechanicalperformance and failure behavior of FRCPCs. Moreover, we validate the pre-sented computational model against experimental data obtained bysubjectingFRCPCstoarangeofthree-pointbendinganddirect tensiletests.Overall,thisresearchoffersaperspectivesforthe de-signandoptimization ofFRCPCs tomature their developmentas functionalandreliableload-bearingbiomaterials.

2. Experimentalstudies

2.1. Three-pointbendingandtensiletests

Poly(vinyl alcohol)(PVA) fiberswere utilizedtoreinforce CPC matrices in order to evaluate how fiber-matrix interfacial prop-erties affect the macro-mechanical properties of fiber-reinforced CPCs. To perform the experimental three-point bending and di-rect tensile test, PVAfiber-reinforced CPC specimens were fabri-catedby firstmixing 98_.5 %

α

-tricalcium phosphate(

α

-TCP, D50 of2.97

μ

m,D90 of6.06

μ

m, volumemeandiameterof 3.51

μ

m) (CAM Bioceramics, Leiden, the Netherlands) with 2.5 wt% PVA fibers (Kuraray Europe GmbH, RFS400_/18mm_, Hattersheim am Main, Germany,200

μ

m indiameter and 18 mm in length) that were manually cut into lengths of8 and 4 mm. Subsequently, a 4wt%NaH2PO4.2H2O(Merck,Darmstadt,Germany)aqueous solu-tion wasaddedat aliquid-to-powder ratio(L/P) of1:2 andfully mixed for 1min until a cementitious paste with randomly dis-persedPVAfiberswasformed [60,61].Theinitialandfinalsetting timesofthesecementformulationwere_{∼ 3}and_∼20min, respec-tively [62,63].Thepastewassubsequentlycastinto polydimethyl-siloxane (PDMS) rectangular molds

(

40× 10× 10 mm

)

, clamped between two glass slides and allowed to set at room tempera-turefor24h.Afterwards,the specimenswere removedfromthe moldsandincubatedinaphosphate-bufferedsalinesolution(PBS) at 37◦ C for 72h to allow the CPC to fully cure. The weight of eachcoatedPVAfiber(l_f=18 mm)was0.77mg,andeach fiber-reinforcedCPCsamplecontained 125mgofPVAfibers. The den-sityofPVAfiberwas∼ 1.3g/cm3 _{andthefibercontentwaskept} constantat2.5wt%forboththethree-pointandtensiletests irre-spective ofthe fiberlength. As a result, forboth the three-point bending and tensile tests, the number of PVA fibers within the each fiber-reinforced CPC sample wascalculated as 370 and 740 forfiberlengths8and4mm,respectively.

To ensure a controlled crack initiation and propagation pro-cess duringthe mechanical tests, a 3_{× 1} mm singleedge notch and two opposing 2.5× 1 mm notches were cut into the speci-mens using a diamond-tipped circular saw blade. These notches prevent an uncontrolled and random failure process and allows thevalidationofthecomputationalmodelparametersaccordingto theexperimental data.Toensurethatall experimental testswere performedunderhydrated conditions,thespecimenswere stored in 37◦ C (PBS) until the moment of testing. Schematic sketches of both three-point bending and tensile test setups, boundary

Fig. 1. Experimental force-displacement curves of (a) the three-point bending test and (b) the tensile test for rectangular FRCPC specimens (40 × 10 × 10 mm ) containing 8 mm PVA fibers ( d f = 200 μm ).

conditions andspecimengeometry are illustrated inFigure SM-1 (seesupplementarymaterial).

Using auniversaltestingmachine(LLOYDmaterialtesting, LS1 series) equippedwith a 1000 N load cell, the three-point bend-ing testandthetensiletestwereperformedatacrossheadspeed of1mm/min.Toperformthethree-pointbendingteststhe spec-imen was placed on two supporting pins with a span length of 30 mmand theloadwasapplied incrementally atthe mid-span. Valuesoffracturetoughness,KIC,werecomputedusing Eq.(1)and (2) based on ASTM C1421− 10 (an established test method for fracture toughness of advanced ceramics at ambient tempera-ture [64])as: KIC=g



Fmaxs010−6 bw3/2



3

(

a/w

)

1/2 2

(

1− a/w

)

3/2



(1) g=g



a w



= 1.99−

(

a/w

)(

1−

(

a/w

))(

2.15− 3.93

(

a/w

)

+2.7

(

a/w

)

2

)

1+2

(

a/w

)

(2) where Fmax is the maximum force and s0 and b represent the support spanandthespecimenthickness, respectively. Thenotch depthisequaltoa, wisthewidthofthetestspecimenandg is

afunctionofratioa/w,asshownin Eq.(2).Thebendingstrength

f_{f l} wascalculatedusing Eq.(3):

ff l=

3s0Fmax

2bw2 (3)

where w denotesthe distancebetweenthe tipof thenotch and thetopofthespecimen.

For tensile tests, the double-notched specimens was first mounted to the clamps of the universal testingmachine. To this end, two plasticT-shapedbars were3-D printedandgluedto ei-ther endof the specimenusing a two-componentPleximon glue (Evonik Röhm GmbH, Darmstadt, Germany). The specimen was subsequentlyplacedintothetestingmachineandsubjectedto ax-ial tensileloadinguntil rupture.Herethetensilestrength fts was calculatedas:

fts=

Fmax

bw (4)

where w is the distance between the tip of the two notches for the tensiletest specimens. The work offracture under bend-ing(WOFb)andtensile(WOFt)loadingweremeasuredbydividing the total area underthe force-displacement curvesby the cross-sectionalsurfacearea [61].

The experimental data were reported as average ± standard deviation and analyzed statistically by means of one-way analy-sis of variance (ANOVA) followed by a Tukey post hoc test. For all tests, force-displacement curves were recorded foratleast 10 samplesperexperimentalgroup(n≥ 10)andadatacollection fre-quency of 16 kHz was applied to record the experimental data points.Thespatialdistributionoffiberswithintestspecimenswas characterizedusingnano-computedtomography(nano-CT,Phoenix NanoTomM,GeneralElectric,Wunstorf,Germany).Nano-CT anal-ysiswasobtainedusingavoxelsize of5.6

μ

m,vocalsize spotof 0.84mm,X-raysourceof 70kV/200

μ

A,andexposuretime of 500mswithouttheapplicationofafilter.

2.2. Resultsanddiscussion

Forthethree-pointbendingtestandthetensiletesttheapplied force was recorded asa function of the vertical displacement of thetop-face of thespecimenat mid-spansection (bending tests) andthe axial displacementof the top-endof thespecimen (ten-siletests).Thecorresponding force-displacementcurvesfor three-pointbendingandtensiletestsofspecimensreinforcedwith8and 4mmPVAfibersareillustratedin Figs.1and 2.

Generally, for both bending and tensile force-displacement curvesof fiber-reinforced CPCs, three main phasescan be distin-guished. In the first phase, the specimen remains in its elastic regime.TheincorporationofPVAfibersenhancethematerial resis-tancetodeformation.Inthesecondphase,byincreasingtheforce, nano- andmicro cracks form near the notch corners where the highest stress concentration occurs. Subsequently, the accumula-tionofthesenano-andmicro-cracksleadstoasinglemacro-crack. Inthethirdphase,tougheningmechanismsbecomeactive and,in the crack wake, PVA fibers partially or fully bridge the crack to stabilizethecrackpropagationprocess.Thebridgingfibersdebond andtransmittheforceacrossacrackandatthefiber-matrix inter-face.Whenfurthercrackpropagationoccurs,fiberpull-outresults in additional frictional sliding resistance, fibrillation of fiber sur-face,andfinally slip-hardeningbehavior. Fiber type,number, dis-tribution,embedded(debonded)lengthandlocationrelativetothe damagezonesignificantlyaffectthepostpeakresponse.Atypical force-displacementcurveobtainedduringthethree-pointbending testofafiber-reinforcedCPCspecimencontaining8mmPVAfibers isdepictedinFigureSM-2(seesupplementarymaterial).

Comparedtothespecimens reinforcedwith8 mmPVAfibers, 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 micro-mechanicalpull-outtest,thecriticalembeddedlength

(

lc

)

ofPVA

Fig. 2. Experimental force-displacement curves of (a) the three-point bending test and (b) the tensile test for rectangular FRCPC specimens (40 × 10 × 10 mm ) containing 4 mm PVA fibers ( d f = 200 μm ).

Fig. 3. Two experimental fiber-reinforced CPC specimens after performing a three-point bending test (left) and a tensile test (right), l e = 8 mm , d f = 200 μm .

Table 1

A summary of experimental results.

Type of test le df KIC ff l fts Ef l Ets WOF b WOF t Number of

mm μmm MPa m 1 /2 _{MPa MPa GPa GPa kJ/ m} 2 _{J/ m} 2 _specimens Three-point bending 8 mm 200 0 . 224 ± 0 . 05 4 . 08 ± 0 . 78 – 1 . 04 ± 0 . 26 – 2 . 35 ± 0 . 56 – 10 Tensile 8 mm 200 – – 1 . 25 ± 0 . 25 – 0 . 16 ± 0 . 08 – 1 . 73 ± 0 . 55 10 Three-point bending 4 mm 200 0 . 186 ± 0 . 02 3 . 34 ± 0 . 48 – 0 . 91 ± 0 . 17 – 0 . 57 ± 0 . 16 – 10 Tensile 4 mm 200 – – 1 . 17 ± 0 . 37 – 0 . 16 ± 0 . 08 – 0 . 44 ± 0 . 06 10

fibers with 200

μ

m diameter embedded in a CPC matrix was 8 mm [51]. Thisimplies that PVA fiberswith 8 mm embedded length arethelongestfiberthatcanbepulledoutfromtheir sur-roundingCPCmatrixwithoutrupture.Thisresultsinhighervalues ofpull-outwork,alongerinterfacialfrictionalresistancephaseand thus more energy dissipation.Consequently, the fiberlength is a prominentparameter that canstrongly affectthebehavior ofthe post-peak regime and the efficiency of the interfacial properties in fiber-reinforced CPCs. Moreover, the presence ofthe fibers af-fectedthe damagewidthandpropagationpattern. Imagesoftwo representative samplesafterexecutionofthebending andtensile tests are shownin Fig. 3.Compared tothe rapid failure induced byanarrowsinglemacro-crackasobservedinfiber-freeCPCtests, we observed a slowercrackgrowth process andawider damage width [53].

A summary of the obtained experimental results is reported in Table1. Inthistable le anddf arethe fiberembeddedlength and fiber diameter, respectively. The fracture toughness, bending strength, tensile strength, bending and tensile modulus are de-notedasKIC, ff l, fts,Ef l andEts respectively.ThevaluesofWOFb andWOFt representtheworkoffracture underbendingand ten-sileloading,respectively.

Fracturetoughness(see Eq.(1))isaquantitativeparameterthat expressesthematerial’sresistancetocrackpropagation.Compared to the previously reported value of fracture toughness for fiber-free CPCmatrixofsimilar chemicalcomposition (0.17MPam1/2_), PVA fibers were able to increase the fracture toughness of fiber-reinforced CPC matrices during the crack initiation and growth

around 9% for 4 mm fibers and 30% for8 mm fibers [53]. This low-foldincreasemightberelatedtosub-optimalfiberdimensions (low fiber aspect ratio for effective reinforcement) which were necessary forthis studyto be able to perform combined experi-mentalandnumericalanalysis.The valuesmeasuredforthe frac-turetoughnessofFRCPCmatricesfitwithintherangereportedfor trabecular bone (0_.1_{− 0.}8 MPam1/2_{), although they are still} be-low valuesreportedforcorticalbone (2− 12MPam1/2_). Parame-terssuch asthe interface affinity, fiberlength,fiber volume frac-tion,fiberalignmentanddistributiondirectlyaffectthe mechani-calpropertiesofFRCPCs [47,65,66].Tailoring theseparameters us-ing computational models asdeveloped herein will contribute to optimizationof themechanical propertiesof fiber-reinforced CPC matrices,whichmightultimatelyleadtobiomaterialswithsimilar propertiesasfoundinbone.

Moreover,incomparisontofiber-freeCPCspecimens,thework offracturevaluesweresignificantly enhanced [47,53].Thismeans that more energy was dissipated during the fracture process of fiber-reinforced specimens, due to the crack-bridging and fiber pull-outprocess.

Thenano-CTimagesofspecimensreinforcedwith8and4mm PVAfibersarepresentedin Fig.4.Theimagesshowvarious side-viewsofdifferentcrosssectionsof specimensafterexecuting the mechanicaltests.Thelightgray,darkgrayanddarkredcolors rep-resenttheCPCmatrix,PVAfibersandfracturesurface,respectively. Thenano- andmicro-poreswere homogeneously distributedover thespecimens.PVAfiberswererandomlydistributedwithrespect tolocationsandorientations.Asdepictedin Fig.4-bandc,thePVA

Fig. 4. Nano-CT images of a CPC matrix reinforced with PVA fibers (a-d) with l e = 8 mm and d f = 200 μm and (e-h) with l e = 4 mm and d f = 200 μm . Light gray, dark gray and dark red colors represent the CPC matrix, PVA fibers and fracture surface, respectively.

Fig. 5. Weibull plots of (a) bending strength and (b) tensile strength.

fibersbridgedthemacro-cracksandstoppedthecrackpropagation process. Inaddition,asshownin Fig.4-aande,some ofthePVA fibers near the damage zone were partially pulled out fromthe matrixfollowedbycontinuation ofcrackpropagation. By compar-ing thefiberlengthof(partially)pulledoutfibersinthenano-CT imagesandexperimentalspecimens,itwasconcludedthatnoneof thePVAfibersrupturedduringthefractureprocessandthe macro-crack mainly grewthrough the CPC matrix.This is inagreement withourmicro-mechanicalpull-outstudywhereacritical embed-dedlengthof8mmwasobservedforPVAfibers [51].Hence,the CPC matrix is considered as a homogeneous material, while the PVAfibersareassumedaslinearelasticmaterialinthenumerical studiesdescribedin Section3.

To study the reliability of the measured material strength, Weibull modulus values were calculated for each experimental group [67,68].TheanalysiswasperformedaccordingtotheASTM standardC1239− 07usingthefollowingequation [69]:

ln



ln



1

(

1− Pf

)



=mln

 σ

σ

0



(5) inwhich Pf = i− 0.5 n (6)

where P_f and m are the failure probability at stress

σ

and the Weibull modulus, respectively. The value

σ

0 is the characteristic stressatwhich63%ofthespecimensfail [70].Thisvaluefor spec-imenscontaining8and4mmembeddedfiberswasmeasuredas 4.6and3.6MPaforbendingtest,and1.5and1.2MPaforthe ten-siletest,respectively.Theparametersnandidenotethetotal num-ber of experimental specimens per experimental group and the specimenrankinascending orderoffailure stresses.The Weibull plotsofthebendingandtensilestrengthofCPCmatrixreinforced with8and4mmPVAfibersarepresentedin Fig.5.

AsummaryofWeibullparametersarelistedin Table2.The cor-relationcoefficientispresentedasR2_.

As shown in Fig. 5, the strength values for both three-point bendingtestandtensiletestareapproximatelyWeibulldistributed. Themeasured Weibullmoduli

(

m

)

forreinforcedCPCs arehigher compared to the Weibull moduli

(

m

)

for fiber-free CPCs (9.7 for bending and 3.6 for tensile tests), which indicates that the reli-abilityof fiber-reinforced sampleswashigher [53].Moreover, the reinforced CPCs containing 4 mm PVA fiber displayed a higher Weibull modulus compared to CPCs reinforced containing 8 mm PVAfibers.Thisphenomenoncanbeexplainedbythefactthatthe larger numberof smaller fibers were moreuniformly distributed throughoutthematrixandcreateamorehomogeneousmaterial.

Table 2

A summary of Weibull parameters of bending and tensile strength.

Type of test fiber embedded length Trend-line formula Weibull modulus ( m ) R2 Three-point bending 8 mm y = 10 . 7 x + 0 . 2 10.7 0.92 Three-point bending 4 mm y = 11 . 1 x + 0 . 2 11.1 0.94

Tensile 8 mm y = 8 . 1 x + 4 . 6 8.1 0.90

Tensile 4 mm y = 10 . 6 x + 0 . 2 10.6 0.94

3. Numericalmodeling

In this section a gradient-enhanced damage model combined with a dimensionally-reduced fiber model is described to char-acterize the failure response ofFRCPC under bendingand tensile loading using finite element modeling. Herein, the main goal is to calibratethenumericalmodelparameters accordingto the ob-tained experimental results andpresent a numerical model that canbeefficientlyusedtofullycharacterizethemechanical perfor-manceofFRCPCs.

3.1. Constitutiverelations

We employ the 3-D gradient-enhanced damage model to nu-merically representthepropertiesofCPC matrix,whilethe mod-eling of PVAfiberswas performedby meansof a dimensionally-reduced fibermodel. Furthermore,a three-phase traction separa-tion law(TSL) is employed to connect the matrix and fiber ele-mentsandrepresentsthefiber-matrixinterfacebehaviour [51].

3.1.1. CPCmatrix

TonumericallymodelthepropertiesoftheCPCmatrix,weused theimplicitgradient-enhanceddamagemodelthatwasdeveloped by Peerlings etal. [54]. Thismodel solves the pathologicalmesh dependence, which is inherently present in standard local dam-agemodels.Herein,forsmalldeformations,therateformof stress-strainrelationfortheCPCmatrixreadsas:

˙

σ

m=

(

1− D

)

Hm:

ε

˙m (7)

where

σ

m and

ε

m arestressandstraintensorandHm andD rep-resent the fourth-order elasticitytensor and damage scalar vari-able, respectively. The damage parameterD variesbetween0 for fullyintactmaterial and1forfullyfailedmaterial. Inadditionto theequilibriumequation,thegoverningfieldequationsincludethe diffusionequation(modifiedHelmholtzequation) [71]:

¯e− c

∇

2_{¯e=e˜ (8)}

where

∇

 is the gradient operator and ¯e and e˜ denote the non-localandlocalequivalentstrain,respectively.Theconstantc isthe materialparameter(gradientparameter)ofthedimensionmaterial lengthscale(lm)square.Thenaturalboundaryconditionisusedas follows:



∇

¯e.n=0 (9)

wherenrepresentstheunitoutwardnormal.Thedamage param-eterDin Eq.(7)isanexplicitfunctionofhistoryparameter

κ

.The initial threshold of historyparameter is

κ

i andits evaluation ac-cordingtotheKuhn-Tuckerrelationsisdefinedas [72,73]:

¯e−

κ

≤ 0,

κ

˙

(

¯e−

κ

)

=0,

κ

˙ ≥ 0 (10) The exponentialdamageevolutionlawisusedtodefinethe dam-agebehavior [74,75]:

D

(

κ

)

=1−

κ

i

κ



(

1−

α)

+

α

exp(−β(κ−κi))



₍₁₁₎

in which

α

and

β

are the materialparameters. The parameter

β

corresponds to therateofdamage growth.In ordertodetermine

theequivalentstrainin Eq.(8),themodifiedvon-Misesdefinition wasemployed [54]: ˜ e= k− 1 2k

(

1− 2

ν

)

I1+



1 2k



k− 1 1− 2

ν

I1

2

+

₍

12k 1+

ν

)

2J2

1/2 (12)

whereI1andJ2 arethefirstandsecondinvariantofthestrain ten-soranddeviatoricstraintensor.

3.1.2. PVAfiber

Accordingtothenano-CTobservations,PVAfiberscanbe con-sidered as linear elastic, while the effect of fiber surface fibril-lation during the pull-out process was assumed negligible and wasthereforeignored [51].Thiscanbealsoexplainedby thefact that thefiberlength doesnotexceed thefibercriticalembedded length (le≤ lc) andtherefore thepull-out process occurswithout fiberrupture.Hence,thestress-strainrelationforthePVAfibersis isotropicallylinearelasticandcanbegenerallydefinedas:

σ

f=Hf:

ε

f (13)

whereHf isthe elasticstiffnesstensorofa PVAfiberand

σ

f and

ε

farethestressandstraintensor.

3.1.3. InterfacebetweenPVAfiberandCPCmatrix

Tointroduce theinterfacialpropertiesbetweenthe PVAfibers and CPC matrix, we employed zero-thickness interface elements with a constitutive traction separation law (TSL). The interfa-cial traction tc along the fiberis a function of the displacement jump:

˙

tc=T[[u˙int]] (14)

wheretc=[ts, tt, tn] andT and[[u˙int]] arethe cohesive tangent matrixandtherateformofdisplacementjumpinlocalcoordinate systemattachedtotheinterfacesurface.Thecohesivetangent ma-trixTcanbewritteninmatrixformas:

T=

_T s 0 0 0 Tt 0 0 0 Tn



(15)

whereTs,TtandTn aretheshearandnormalstiffnessesofthe in-terface.Herein,weonlyconsiderfiberslipalongthefiberaxis di-rection,therebythevaluesofTt andTn,stiffnessesoftheinterface inthedirectionnormalandperpendiculartotheinterfacesurface andfiberaxis,aresetartificiallyhightoprevent interface separa-tioninthosedirections.Thedisplacementjumpdenotetheslip be-tweenfiberandmatrixsintangentialdirection.Thepreciseform oftheinterfacetangentialstiffnessTs=

∂

ts/

∂

sdependsontheTSL used. In our previous work [51], we performedmicromechanical pull-out experimentsto investigate the affinity betweenthe PVA fibersandthe CPC matrix.Asa result, we proposed theuse ofa three-phaseTSL tonumericallymodelthe completepull-out pro-cessi.e. theelastic, debondingandfrictional pull-outphases. Ac-cording to the three-phase TSL formulation, the bond stress-slip

Fig. 6. Problem domain and boundary .

relationreadsas:

ts

(

s

)

=

⎧

⎪

⎨

⎪

⎩

Gs 0≤ s≤ s0

τ

maxexp



s0−s df γ0



s0≤ s≤ s1 df



γ

2

(

s− s1

)

2+

γ

1

(

s− s1

)



+

τ

maxexp



s0−s1 df γ0



s>s1 (16) in which G represents the corresponding relative bond modulus and

γ

0,

γ

1, and

γ

2 denote the parameters controlling the de-scending andascending branchesinthe debondingandfrictional stages.The startingpointofthe debondingphaseandthesliding phase are presented ass0 ands1, respectively(see [51]for more details).

3.2. Weakformofgoverningequations

Theweakformsofequilibriumanddiffusionequations govern-ingthedeformationandfailureprocessesofafiber-reinforced cal-ciumphosphatecementspecimenare presentedbelow.The prin-ciple ofvirtualworkinadomain



withtheboundaries

ofthe domaincanbedescribedasfollows(see Fig.6):

δ

Wm







matrix +

δ

Wf



fiber +

δ

Wintf







interface







internalforces =

δ

Wext







externalforces (17) with

δ

Wm₌  m

∇

s

_δ

_u m:

σ

md



m (18)

δ

Wf₌  f

∇

s

_δ

_u f:

σ

fd



f (19)

δ

Wintf₌  int

δ

[[uint]].tcd

int (20)

δ

Wext=  

δ

u.fbd



+  t

δ

u.¯td

t (21)

inwhich,

δ

uisthevirtualdisplacementappliedovertheDirichlet boundary

u

andthetraction¯timposedovertheNeumann bound-ary

t (

=

u∪

t).Thebodyforcevectorispresentedasfb and

∇

s _{is the symmetric gradient operator. The virtual displacement} vectors over the body domain



=

m

∪

f

formatrix and fiber are denotedas

δ

um and

δ

uf,respectively.Thedisplacementjump andinterfacialtractionsaredefinedover

int

whichcorrespondsto

theinterfacesurfaceofallembeddedfibers [55,76].Theweakform ofthediffusionequation(see Eq.(8))canbeformulatedas: 

m



δ

¯e¯e+c

(

∇δ

¯e

)

T

∇

¯e



d



m= 

m

δ

¯ee˜d



m (22) Eqs.(17)and (22)formthegoverningsystemofcoupledsetof par-tialdifferentialequations.Inthefollowing,thegeneralformofthe globalsystemofequationsforthe finiteelement implementation isdescribed.

3.3. Spatialdiscretization

Following the Galerkin approach, the displacement vector for matrix um andfiber uf and the nonlocal equivalentstrain ¯e can bediscretizedbymeansofshapefunctionsas:

um=Nmu

˜m uf=Nfu˜f ¯e=Ne˜¯e (23)

where matrix Nm,Nf and rowvector Ne contain the shape func-tions for the displacement fields and the non-local equivalent strains. Moreover, u

˜m, u˜f and ˜¯e represent the nodal degrees of freedomofthematrixandfiberdisplacementcomponentsandthe nonlocalequivalentstrains,respectively.Foreachinterfaceelement thedisplacementjumpcanbedeterminedas:

[[u˙int]]=ANintu˙gint (24)

where A and ug_int=[u

˜f,1,u˜f,2,u˜m,1,u˜m,2,...,u˜m,n]are the rotation matrixfromtheglobaltothelocalcoordinatesystemandthe dis-placementvectorforfiberandmatrixintheglobalcoordinate sys-tem, respectively. The nodal displacements of a one-dimensional fiberelementaredefinedasu

˜f,1 andu˜f,2 (node1and2)andu˜m.n correspondstothedisplacementatnodenoftheparentmatrix el-ement.Moreover,Nint=[Nint,f,Nint,m]denotestheshapefunctions oftheinterfaceelements,inwhich:

Nint,f=

⎡

⎢

⎣

Nf,1 0 0 Nf,2 0 0 0 Nf,1 0 0 Nf,2 0 0 0 Nf,1 0 0 Nf,2

⎤

⎥

⎦

(25) Nint,m=

⎡

⎢

⎣

Nm,1 0 0 ... Nm,n 0 0 0 Nm,1 0 ... 0 Nm,n 0 0 0 Nm,1 ... 0 0 Nm,n

⎤

⎥

⎦

(26) whereNf,1 andNf,2 aretheshape functionsofa one-dimensional fiber element (node 1 and 2) and Nm,n is the shape function at node n of the parent matrix element [55]. Differentiation of

Eq.(23)leadstothestrain componentsandgradientofnon-local equivalentstrain:

ε

m=Bmu

˜m

ε

f=Bfu˜f

∇

¯e=Be˜¯e (27) whereBm,BfandBe containthederivativesoftheshapefunctions

Nm,NfandNe,respectively.

3.4. Linearizationandtheincremental-iterativesolutionprocedure

The next step consists of the linearization of the governing equationsin order to constructa consistent incremental-iterative Newton-Raphson solution procedure. In this regard, linearization atiterationi+1atthenodallevelleadsto:

u ˜ i+1 m =u ˜ i m+

u ˜ i+1 m u ˜ i+1 f =u_˜ i f+

u_˜ i+1 f ¯e_˜ i+1_=¯e ˜ i +

¯e ˜ i+1 ₍₂₈₎

andatintegrationpointlevelwehave:

ε

i+1 m =

ε

im+

ε

mi+1

ε

if+1=

ε

i f+

ε

if+1

σ

i+1 m =

σ

mi +

σ

mi+1

σ

fi+1=

σ

i f+

σ

fi+1 (29) Di+1_=Di₊

_Di+1 _e˜i+1_=e˜i₊

_e˜i+1 inwhich,

ε

i+1 m =Bm

u ˜ i+1 m

ε

if+1=Bf

u ˜ i+1 f

σ

i+1 m =

(

1− Di

)

HmBm

˜u i+1 m − Hm

ε

im

Di+1

σ

fi+1=HfBf

˜u i+1 f

Di+1=



∂

D

∂κ



i



∂κ

∂

¯e



i Ne

¯e ˜ i+1

_e˜i+1₌



∂

e˜

∂ε

m



i Bm

u ˜ i+1 m (30)

where

representstheiterativeincrementbetweeniterationiand

i+1.In Eq.30,therelationbetweenthehistoryparameter

κ

and ¯e(

∂

κ

_/

∂

¯e)isequalto1when ¯e_>

κ

0and0otherwise.Herein,

κ

0 de-notestheconvergedvalueofthehistoryparameterintheprevious increment.Thesetofgoverning Eqs.(17)and (22)atiterationi+1 canbeexpressedas:

f_inti+1=fext− finti (31)

finti+1,e=fext,e− finti ,e (32) wherefextandfint denotethediscrete balanceofexternaland in-ternal nodal forces,respectively. Similar expressions are used for thefext,e,fint,einthediffusionequationforclarity(fext,e=0). Sub-stitution of expressions 28–30into Eqs. 31 and (32)and follow-ingthestandardlinearizationproceduregives(moredetailscanbe foundin [54,55,77]):  m BT_m



1− Di



_H mBm

u ˜ i+1 m d



m−  m



∂

D

∂κ



i



∂κ

∂

¯e



i × BT mHm

ε

miNe

¯e ˜ i+1_d

_

m+  f BT_f

σ

i f

u_˜ i+1 f d



f (33) +  int

ATNTinttci

[[uint]]i+1d

int=fext− finti

 



NTeNe +cBTeBe



¯e ˜ i+1_d

_

−  m NTe



∂

e˜

∂ε

m

T

i Bm

u ˜ i+1 m d



m=−finti ,e (34)

The last step to construct the global systemof equations,is in-troducing the stiffnessmatrices to rephrase the set ofgoverning equationsintoacompactmatrixform.

3.4.1. Matrix

ConsideringthetermsrelatedtotheCPCmatrixinEqs.(33)and

(34),thefollowingstiffnessmatricescanbeintroduced [72]:

Kmm=  m BT_m



1− Di



_H mBmd



m (35) Kme=−  m



∂

D

∂κ



i



∂κ

∂

¯e



i BT_mHm

ε

miNed



m (36) Kem=−  m NTe



∂

e˜

∂ε

m

T

i Bmd



m (37) Kee=  



NTeNe +cBTeBe



d



(38) 3.4.2. Fiber

Accordingtothedimensionally-reducedfibermodel [55,56],all PVA fibers were considered asone dimensional objects that can transfertheloadonlyinaxial direction.Fibercouplingis ignored forsimplicity.Inaddition,since thefibersare discretizedwith 1-D elements, overlap of fibers is not considered. This idealization offibershelpstodrasticallyreducethecomputationalcostduring the numerical simulations compared to the conformal finite ele-mentapproach,however,stillsixextradegreeoffreedomperfiber elementareadded.Thevalidityofthedimensionally-reducedfiber modelcompared to a conforming3-D fibermesh isdiscussed in detail in [55]. Using the dimensionally-reduced fiber model, the PVAfiber discretizationis independent fromtheCPC matrix dis-cretizationand,moreover,therearenoconstraintsonthenumber ofPVAfibersthat intersect aCPC matrixelement.This simplifies themeshdiscretizationprocess.ThePVAfibersarediscretized us-ingtypicaltwo-noded3-D trusselementsinwhicheachnodehas three globaldegrees offreedom. To simplifythe calculation pro-cessofthestiffnessmatrixformatrixelementsthatarecrossedby fibers,thematrixcontributioniscalculatedoverthetotal volume oftheelement.ThereforeaneffectiveYoung’smodulus isusedto generatethefiberstiffnessmatrixandcanceloutthealready com-putedmatrixcontributioninthefiberregion.The stiffnessmatrix foreachfiberelementisgivenby [78]:

(39)

whereE_f,Em,Af andLef denotethefiberandmatrixYoung’s mod-uli,fibercrosssectionarea andfiberelement length,respectively. Usingthenodecoordinates,thecosinescanbedefinedas:

Cx=x2− x1 Le f Cy= y2− y1 Le f Cz=z2− z1 Le f (40) inwhich theglobalcoordinate offiberelement endsare defined as(x1, y1, z1)and(x2, y2, z2).

3.4.3. Fiber-matrixinterface

Thecohesive tangentstiffnessmatrixforagiveninterface ele-mentbetweenPVAfiberandCPC matrix(see Eqs.(24)and (33))

canbeformulatedas [55,77]: Kint=



Kint,mm −Kint,mf −Kint,fm Kint,ff



(41) with Kint,mm=  int

NT_int,mATTANint,md

int (42)

Kint,mf= 

int

NTint,mATTANint,fd

int (43)

Kint,fm= 

int

NTint,fATTANint,md

int (44)

Kint,ff= 

int

NTint,fATTANint,fd

int (45) and

d

int=CfdLef (46)

Considering fibers as one-dimensional objects, for each sub-matrix in Eq.(41),the integrationwascomputedover thelength ofinterface elementLe

f (equaltofiberelementlength)and multi-pliedbyC_fthatdenotesthefibercross-sectioncircumference.

3.5. Globalsystemofequations

Byassemblingallsub-matricesintroducedinprevioussections, theglobalstructureofcoupledsystemofequations(Eqs.(33)and

(34))canbewrittenas:

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎣

K_mm K1_mf K_mf2 ... Kn_mf Kme K1fm K 1 ff 0 ... 0 0 K2fm 0 K2ff ... 0 0 . . . ... ... ... ... ... Knfm 0 0 ... Knff 0 Kem 0 0 ... 0 Kee

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎦

⎡

⎢

⎢

⎢

⎢

⎢

⎣

u ˜m

u ˜ 1 f

u ˜ 2 f . . .

u ˜ n f

e ˜

⎤

⎥

⎥

⎥

⎥

⎥

⎦

=

⎡

⎢

⎢

⎢

⎢

⎢

⎣

fext 0 0 . . . 0 0

⎤

⎥

⎥

⎥

⎥

⎥

⎦

−

⎡

⎢

⎢

⎢

⎢

⎢

⎢

⎣

fint,m f1int,f f2int,f . . . fn_int,f fint,e

⎤

⎥

⎥

⎥

⎥

⎥

⎥

⎦

(47)

inwhich,Kmm=Kmm+Kint,mm,Kmf=Kint,mf,Kfm=Kint,mf,Kff=

K_int,ff+K_fand: fext=  N T mfbd



+  t NT_m¯td

t (48) fint,m=  m BT_m

σ

i md



m+Cf  Le f NT_int,mATTANint,mu ˜mdL e f −Cf  Le f NTint,mATTANint,fu ˜fdL e f (49) fint,f=−Cf  Le f NTint,fATTANint,mu ˜mdL e f +Cf  Le f NTint,fA T TANint,fu ˜fdL e f+  f BTf

σ

fid



f (50) fint,e=  



NT_eNe +cBTeBe



¯e ˜ i_d

_

₋ N T ee˜id



(51)

where the contribution of interface stiffness is shared between fiber andmatrix degree of freedom. Eq. (47) is solved using an implicitNewton-Raphsonscheme.Ineach iteration,thesystemof

equationsaround theapproximatesolution calculatedinthe pre-viousiterationislinearizedandsolvedtocomputeanew approx-imation of the solution.This procedure is repeated until conver-genceofthe solutionis reached. Herein,theiteration superscript

i+1andirelatedtotheimplicitNewton-Raphsonsolverare omit-tedforclarity.Theproposednumericalmodelforsolvingthis non-linear system of governing equations is coded in the research-orientedC++programmingtoolkitnamedJive [79].Inthenext sec-tion, the results of the numerical simulations are presented and themodelisvalidatedagainsttheexperimentalresults.

3.6. Resultsanddiscussion

A3-D FEmodel wasdeveloped tonumerically studythe fail-ure process offiber-reinforced calciumphosphate cements under three-pointbending testandtensiletest.A3-D modelrepresents the fiber distribution in the CPC discretization more realistically despite the fact that it is computationally more costly compared toatwo-dimensional(2-D)numericalframework.

Inallnumericalsimulations theCPCmatrixwasdiscretized by tetrahedralcontinuum elementsandtrusselementswereused to discretize the PVAfibers. The PVA fibers remain elastic and will notbreakduringthepull-outprocess.Moreover,thezero-thickness interfaceelementswithathree-phaseTSLwere incorporatedinto themodeltorepresenttheinterfacialpropertiesbetweenthePVA fibersandCPCmatrix.Theproblemboundaryconditionswere en-forced asillustrated in Figure SM-1. For the three-point bending test, the specimen is subjected to a prescribed displacement at mid-span,theleftsupportingpinisfixedandtherightpinis con-strainedin y-andz-directions.For thetensiletest, thespecimen issubjectedtoaxialtensileloadingwhilethe bottomedgeofthe specimenisfixed.

Tomimicthe effectsofPVAfibers beingrandomlydistributed duringthe experimental tests, the PVAfibers with uniformfiber densitydistributionarearbitrarydistributedintheCPCmatrix dis-cretization.Foreachnumericalsimulation,tworandomfiber distri-butionwereconsidered.Furthermore,toachieveanaccurate mea-surementof thestress profile over thevicinity of notcheswhere thecrackinitiationanddamagepropagationareexpected,we dis-cretizedthisareawithafinermeshdiscretizationcomparedtothe restofthestructure.

Themeshdiscretization forboth three-pointbendingand ten-siletestsoffiber-reinforced calciumphosphatecementspecimens withtwo random PVAfiber distributions andtwo differentfiber lengths(le=8and4mm)anddiameterdf=200

μ

maredepicted in Figs. 7 and 8, respectively, where h represents the element sizeintherefinedregion.Inthesefigures,forclarity,justthe out-lineandhalfofthematrixare illustratedforthefirstandsecond fibersdistribution,respectively.Moreover,toshowthedifferences, both fiber distributions, black andred lines, are depicted on top ofeach other(2-Dview)forbothbendingandtensilesimulations (seepartseandfin Figs.7and 8).

Thenumericalforce-displacementcurvesforthree-point bend-ingandtensiletestsoffiber-reinforcedCPCwithtworandomPVA fiberdistributionsandtwodifferentfiberembeddedlengths(le=8 and4mm) together withthecorresponding experimental results are plottedin Figs. 9 and 10respectively. We selectedthe three-point bending test forcalibration purposes due to themore sta-ble and controlled crack propagation duringfailure as compared totensiletesting [53].Inthisregardthenumericalmodel parame-terswere tuned in orderto interpolatethe average ofmaximum 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 mechanicalresponseofthefiber-reinforced CPCunderthetensile tests. A B-spline interpolation technique wasused to find a

rep-Fig. 7. Mesh discretization and fiber distributions of the three-point bending test (a, b) and the tensile test (c, d). Both random fiber distributions, black and red lines, are shown in (e) for the three-point bending test and (f) for the tensile test (notched rectangular FRCPC specimens 40 × 10 × 10 mm with PVA fibers l f = 8 mm , d f = 200 μm and h = 0 . 05 mm ).

Fig. 8. Mesh discretization and fiber distributions of the three-point bending test (a, b) and the tensile test (c, d). Both random fiber distributions, black and red lines, are shown in (e) for the three-point bending test and (f) for the tensile test (notched rectangular FRCPC specimens 40 × 10 × 10 mm with PVA fibers l f = 4 mm , d f = 200 μm and h = 0 . 05 mm ).

resentativeaveragecurve foreach setofexperiments(dashedred lines in Figs. 9 and 10) andthisiscompared withthe numerical model.

The numerical simulations were performed with damage threshold

κ

_i=0_.0035_,gradientparameterc₌0_.02mm2_,Poisson’s ratioformatrixandfibers

ν

=0.2andk=18inthemodified von-Misesdefinition(see Eq.(12)) [53].Thedamageevolutionvariables

α

and

β

in Eq.(11) areset as0.9and80,respectively.Moreover, the set of parameters for the interfacial traction separation law

were employedaccordingto thetuned parametersobtainedfrom micro-mechanicalpull-outtestspresentedin [51].

Thenumericalresultsfitwithintheexperimentalenvelope.This evidences that the proposed numericalmodel can reliably repre-sent the bending and tensile responses of fiber-reinforced CPCs. The numerical results of both fiber distributions for each three-pointbending andtensilesimulations arerelativelysimilar. Small differences,mainlyinthepostpeakregime,arerelatedtothe ran-domfiberlocationsandthereforetoadifferentleveloffiber

pull-Fig. 9. Force-displacement curves of (a) the three-point bending test and (b) the tensile test on rectangular FRCPC specimens (40 × 10 × 10 mm ) with 8 mm PVA fibers.

Fig. 10. Force-displacement curves of (a) the three-point bending test and (b) the tensile test on rectangular FRCPC specimens (40 × 10 × 10 mm ) with 4 mm PVA fibers.

Fig. 11. The fully damaged profile for (a) three-point bending test and (b) tensile test (notched rectangular FRCPC specimens 40 × 10 × 10 mm with PVA fibers l e = 8 mm and d f = 200 μm ).

Fig. 12. The fully damaged profile for (a) three-point bending test and (b) tensile test (notched rectangular FRCPC specimens 40 × 10 × 10 mm with PVA fibers l e = 4 mm and d f = 200 μm ).

out that occurs within the damage zone. The damage profile for the three-point bending test and the tensile testsof CPC matrix reinforcedbyPVAfiberswithdiametersof200

μ

mandembedded lengthsof8and4mmarepresentedin Figs.11and 12.The dam-ageprofilesofbothrandom fiberdistributions forthethree-point bendingtestandthetensiletestswerealmostidentical.Therefore,

thedamageprofilesofthefirstrandomfiberdistributionare pre-sentedasrepresentativeillustrations.

For both three-point bending and tensile tests of fiber-reinforced CPCs, the damage profiles are wider compared to the fiber-free CPCs [53]. This was also observed experimentally by comparingthefracturesurfaceofspecimensafterexecutionofthe

tests. For fiber-free CPCs, a narrow and localized single macro-crack with a smooth fracture surface was detected, whereas for fiber-reinforced CPCs a relativelywider and expandedcrack with atortuousfracturesurfacewasobserved(see also Fig.3).In addi-tion,thedamageprofileoftheCPCmatrixreinforcedbyPVAfibers withanembeddedlengthof8mmiswidercomparedtothose re-inforcedwithafiberembeddedlength of4mm (see Figs.11and

12).Duetothecrack-bridgingprocess,thecrackcanzipalongthe fibers andmore fracture energy can be dissipated as a result of fiberpull-out.Thiscanleadtoawiderdamagezone,largervalues ofworkoffracture,andamorestablefailureprocess.

Another parameter that can affect the width of the damage band isthe materiallength scale,whichis relatedtotheaverage sizeofcoalescedporesandmaterialmicro-structureimperfections anditseffectwasintroducedbygradientparameterc inthe diffu-sion equation tothesystemof governingequations.Theeffect of thisparameterforCPCmatrixisstudiedinourpreviouswork [53]. Infact,thegradientparameterc in Eq.(8)adjuststhewidthofthe damageband.Highervaluesofthegradientparametercorrespond to a larger size of micro-structural voids andimperfections that leadtoawiderlocalizationzone.Furthermore,thematerial param-eter

β

inthedamageevolutionlaw(Eq.(11))determinesthecrack propagation rateandlower values of

β

results ina slower crack propagationandconsequentlyamoreductilefailureresponse.The improved mechanical properties of reinforcedcomposites are at-tributedtothehighernumberofinterfaces,therebyincreasingthe numberofchannelsforcrackpropagation.Inaddition,the forma-tion of additional microstructuraldefects produced by immature fiber-matrix interfaces facilitates crack deflection which leads to enhanced fracture energy dissipationand preventsthecomposite fromcatastrophicfailure.Hence,tomimicthewiderdamageband andmorestabledamagegrowththatwasexperimentallyobserved forfiber-reinforcedCPC,weselectedalargermateriallengthscale

c andlowermaterialparameter

β

foroursimulationscomparedto ourpreviousstudyonfiber-freeCPCs [53].

Thegoodagreementbetweenthenumericalresultsandthe ex-perimental envelope evidences that the proposed combined nu-merical modelcanbereliably appliedtofurtherinvestigationson FRCPCs andreplace the tedious experimental trial–and–error de-sign procedures that are commonly used in biomaterials science andengineering.

4. Conclusions

In this work, the mechanical response of a CPC matrix rein-forcedbyPVAfiberssubjectedtobendingandtensileloadingwas investigated both experimentally and numerically. To develop a comprehensivenumericalmodel,weproposed theuseof3-D gra-dientenhanceddamagemodelcombinedwithadimensionally re-duced fiber modelto numerically characterize the failure behav-ior of a fiber-reinforced CPC matrix.The PVAfibers are idealized asone-dimensional objectswhichdrasticallyreducedthe compu-tational cost during the numerical simulations compared to the conformalfiniteelementapproach.Moreover,anadvancedtraction constitutive law is employed to represent the fiber-matrix inter-facial properties.We tuned the parameters of the combined nu-merical modelaccordingtothe averageofexperimental dataand showedthat theproposedmodelisableto predictthefailure re-sponseofafiber-reinforcedCPCmatrixunderbendingandtensile loadingwithgoodaccuracy.WeshowedthatthePVAfiber embed-dedlengthisa keyparameterthatcanaffecttheamountof frac-ture energydissipation,damage zone andstabilityoffailure pro-cess. Wearguethat theproposedcombinednumericalmodelcan beemployedwithagoodapproximationforfurtherstudiesonthe developmentandoptimizationoffiber-reinforcedCPCs.

DeclarationofCompetingInterest

Theauthorsdeclarethattheyhavenoknowncompeting finan-cialinterestsorpersonalrelationshipsthatcouldhaveappearedto influencetheworkreportedinthispaper.

Acknowledgments

Theresearchleadingtotheseresultswasfundedbythe Nether-lands Organization for Scientific Research (NWO, VIDI project 13455).

Supplementarymaterial

Supplementary material associated with this article can be found,intheonlineversion,atdoi:10.1016/j.actbio.2020.10.014.

References

[1] M. Bohner , L. Galea , N. Doebelin , Calcium phosphate bone graft substitutes: Failures and hopes, Journal of the european ceramic society 32 (11) (2012) 2663–2671 .

[2] J. Zhang , W. Liu , V. Schnitzler , F. Tancret , J.-M. Bouler , Calcium phosphate ce- ments for bone substitution: chemistry, handling and mechanical properties, Acta Biomaterialia 10 (3) (2014) 1035–1049 .

[3] M.-P. Ginebra , M. Espanol , E.B. Montufar , R.A. Perez , G. Mestres , New process- ing approaches in calcium phosphate cements and their applications in regen- erative medicine, Acta Biomaterialia 6 (8) (2010) 2863–2873 .

[4] M. Bohner , Calcium orthophosphates in medicine: from ceramics to calcium phosphate cements, Injury 31 (20 0 0) D37–D47 .

[5] W. Habraken , P. Habibovic , M. Epple , M. Bohner , Calcium phosphates in biomedical applications: materials for the future? Materials Today 19 (2) (2016) 69–87 .

[6] W. Liu , J. Zhang , G. Rethore , K. Khairoun , P. Pilet , F. Tancret , J.-M. Bouler , P. Weiss , A novel injectable, cohesive and toughened si-hpmc (silanized-hy- droxypropyl methylcellulose) composite calcium phosphate cement for bone substitution, Acta biomaterialia 10 (7) (2014) 3335–3345 .

[7] M. Epple , K. Ganesan , R. Heumann , J. Klesing , A. Kovtun , S. Neumann , V. Sokolova , Application of calcium phosphate nanoparticles in biomedicine, Journal of Materials Chemistry 20 (1) (2010) 18–23 .

[8] L.L. Hench , J. Wilson , Introduction, in: An introduction to bioceramics, World Scientific, 1993, pp. 1–24 .

[9] E. Sanzana , M. Navarro , F. Macule , S. Suso , J. Planell , M. Ginebra , Of the in vivo behavior of calcium phosphate cements and glasses as bone substitutes, Acta Biomaterialia 4 (6) (2008) 1924–1933 .

[10] I. Ajaxon , A. Acciaioli , G. Lionello , M.-P. Ginebra , C. Öhman-Mägi , M. Baleani , C. Persson , Elastic properties and strain-to-crack-initiation of calcium phos- phate bone cements: Revelations of a high-resolution measurement technique, Journal of the mechanical behavior of biomedical materials 74 (2017) 428–437 .

[11] S.S. Bajammal , M. Zlowodzki , A. Lelwica , P. Tornetta III , T.A. Einhorn , R. Buckley , R. Leighton , T.A. Russell , S. Larsson , M. Bhandari , The use of calcium phosphate bone cement in fracture treatment: a meta-analysis of randomized trials, JBJS 90 (6) (2008) 1186–1196 .

[12] S. Larsson , T.W. Bauer , Use of injectable calcium phosphate cement for fracture fixation: a review, Clinical Orthopaedics and Related Research® 395 (2002) 23–32 .

[13] C.D. Friedman , P.D. Costantino , S. Takagi , L.C. Chow , Bonesourcehydroxyapatite cement: a novel biomaterial for craniofacial skeletal tissue engineering and re- construction, Journal of biomedical materials research 43 (4) (1998) 428–432 .

[14] B.R. Constantz , I.C. Ison , M.T. Fulmer , R.D. Poser , S.T. Smith , M. VanWagoner , J. Ross , S.A. Goldstein , J.B. Jupiter , D.I. Rosenthal , Skeletal repair by in situ for- mation of the mineral phase of bone, Science 267 (5205) (1995) 1796–1799 .

[15] J.S. Al-Sanabani , A .A . Madfa , F.A . Al-Sanabani , Application of calcium phosphate materials in dentistry, International journal of biomaterials 2013 (2013) .

[16] A. Uchida , N. Araki , Y. Shinto , H. Yoshikawa , E. Kurisaki , K. Ono , The use of calcium hydroxyapatite ceramic in bone tumour surgery, The Journal of bone and joint surgery. British volume 72 (2) (1990) 298–302 .

[17] W. Liu , J. Zhang , P. Weiss , F. Tancret , J.-M. Bouler , The influence of differ- ent cellulose ethers on both the handling and mechanical properties of cal- cium phosphate cements for bone substitution, Acta biomaterialia 9 (3) (2013) 5740–5750 .

[18] M. Ikenaga , P. Hardouin , J. Lemaître , H. Andrianjatovo , B. Flautre , Biomechani- cal characterization of a biodegradable calcium phosphate hydraulic cement: a comparison with porous biphasic calcium phosphate ceramics, Journal of Biomedical Materials Research: An Official Journal of The Society for Bioma- terials, The Japanese Society for Biomaterials, and the Australian Society for Biomaterials 40 (1) (1998) 139–144 .

[19] A.J. Ambard , L. Mueninghoff, Calcium phosphate cement: review of mechanical and biological properties, Journal of Prosthodontics 15 (5) (2006) 321–328 .