Multiscale modeling of the effect of sub-ply voids on the failure of composite materials

Delft University of Technology

Multiscale modeling of the effect of sub-ply voids on the failure of composite materials

Turteltaub, Sergio; de Jong, Gijs

DOI

10.1016/j.ijsolstr.2019.01.031

Publication date

2019

Document Version

Final published version

Published in

International Journal of Solids and Structures

Citation (APA)

Turteltaub, S., & de Jong, G. (2019). Multiscale modeling of the effect of sub-ply voids on the failure of

composite materials. International Journal of Solids and Structures, 165, 63-74.

https://doi.org/10.1016/j.ijsolstr.2019.01.031

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.

‘You share, we take care!’ – Taverne project

https://www.openaccess.nl/en/you-share-we-take-care

Otherwise as indicated in the copyright section: the publisher

is the copyright holder of this work and the author uses the

Dutch legislation to make this work public.

International Journal of Solids and Structures 165 (2019) 63–74

ContentslistsavailableatScienceDirect

International

Journal

of

Solids

and

Structures

journalhomepage:www.elsevier.com/locate/ijsolstr

Multiscale

modeling

of

the

effect

of

sub-ply

voids

on

the

failure

of

composite

materials

Sergio

Turteltaub

∗

,

Gijs

de

Jong

Delft University of Technology, Faculty of Aerospace Engineering, Kluyverweg 1, 2629 HS Delft, the Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history:

Received 21 November 2018 Revised 11 January 2019 Available online 31 January 2019

Keywords:

Multiscale fracture Composite materials Voids

Representative surface element Cohesive elements

a

b

s

t

r

a

c

t

Amultiscalefracturemodelisdevelopedtostudytheinfluenceofdefectsappearingatamicroscalein afiber-reinforcedcompositelaminate.Themodelestablishesalinkbetweenthegeometrical characteris-ticsofsub-plyimperfectionsthatmaybecreatedduringmanufacturingandtheoverallfracturestrength andfractureenergyofthecomposite.Inparticular,arecently-developedmultiscaletheoryisexpanded toaccountformicrovoidsinsidethematrixandgapsbetweenclosely-spacedfibersthatpreventfilling. Thesedefectsareexplicitlyincorporatedinfiniteelementsimulationstostudytheirinfluenceonthe on-setandpropagationofcracksatthesub-plylevel.Toconnectthesemicrocrackstotheeffectivefracture behaviorataply-level,acomputationalhomogenizationtechniqueisappliedtoextractthe energetically-equivalentmacroscopicfractureproperties.Throughaparametricanalysisofconfigurations,theinfluence ofthevoidcontent(porosity), voidtypeand voidshapeontheeffectivefracturestrengthandthe ef-fectivefractureenergyofacompositearequantified.Resultsshowthattheporosityisthemain param-eterinfluencingfracturepropertieswhiletheshapeofthedefectsand theirtype(matrixorinterfiber) onlyplayasecondaryrole.Furthermore,theinfluenceofvoidsonthefracturepropertiesappearstobe stronglydependentontheloadingconditions.Inparticular,fortherangeofporosityanalyzed(upto8%), theinfluenceofvoidsinmodeIonthetransversefracturestrengthisnotsignificantbutthetransverse fractureenergydecreasesapproximatelylinearlydowntoabout50%ofitsoriginalvalue.Incontrast,in modeII,thetransversefracturestrengthissignificantlyaffectedwithincreasingporosity.Furthermore, thetransversefractureenergydependsnonlinearlyontheporosityandthereductionisrelativelymore pronouncedthanformodeI.

© 2019ElsevierLtd.Allrightsreserved.

1. Introduction

Defectsthatappearduringmanufacturingofcomposite materi-alsmayaffecttheir effective(macroscopic)mechanicalproperties, in particular their fracture properties such as fracture strength andfractureenergyduringquasi-static,fatigueandimpactloading (Bowles and Frimpong, 1992; Jeong, 1997; Summerscales, 1998; Costaetal.,2001;Hamidietal.,2004;Scottetal.,2014;Maragoni etal.,2017).Thesedefectscanrangefromrelativelylargedryspots spanning a ply thickness to matrix microvoids smaller than the diameter of the fibers. Experimentally quantifying the influence of defectson theeffective mechanicalpropertiesis a challenging tasksince,bytheirownnature,itisdifficulttogeneratedefectsin acontrolled fashion(i.e.,changingdefectswhilekeepingall other parameters fixed Costa et al. (2001)). Due to this large experi-mentalvariability,itisusually notpossibletouniquelyascribean

∗ _{Corresponding author.}

E-mail address: S.R.Turteltaub@tudelft.nl (S. Turteltaub).

effectto specificfeatures in defects.However, recently-developed multiscale numerical techniques can assist in providing insight into the influence of defects through computational simulations. Thesetechniqueshavetheadvantage ofbeingabletocontrol the typesofdefectsincludedinthesimulations and, thus,beingable to establish a direct link between specific types of defects and theirinfluenceontheobservable(macroscopic)properties.

The multiscale approach couples the macro and micro re-sponsesusingahomogenizationtechnique,whichreliesonthe be-havior of the fine scale to determine the response of the larger scale.This conceptalso enablesthe developmentofa failure cri-terion based on the micromechanical behavior of the composite, both in terms of initiation and evolution, and has been exten-sively used for composite materials (see, e.g., Xia et al., 2001; Hettich et al., 2008; Canal et al., 2009; Arteiro et al., 2015). In the present work, attention is given to two distinct types of microscopic defects, namely (i) microvoids inside the matrix and (ii) interfiber voids, which are gaps between closely-spaced fibersthat appear when theresin is unable to flow betweenthe fibers.Thesetypesofdefectshavebeenstudiednumericallyusing

https://doi.org/10.1016/j.ijsolstr.2019.01.031

micromechanicalmodelsinVajarietal.(2013,2014);Vajari(2015);

Maragoni et al. (2016) where it was found that the microscopic voidshavealargeinfluenceonboththemacroscopicstrengthand propagation of damage. The present study investigates a similar situationtoprovidefurtherinsightontheinfluenceofthe micro-scopic characteristics of voids on the fracture properties. In par-ticular,two novelaspects areincludedhere:Firstly,fromthe the-oretical point of view, a recently-developed multiscale technique (Turteltaub et al., 2018) is expanded to account for sub-ply (mi-croscopic)defects and a newaveraging techniqueis proposed to satisfythe scale transitionrequirements(crack-basedHill-Mandel condition).Secondly,asystematicparametricanalysisiscarriedout untilcompletefailureofthe(numerical)specimensand,througha post-processingtechnique,thedetailedresultsarecondensedinto effectivematerialpropertiesforthatmicrostructure.Through para-metricanalyses it ispossible to establish a relation between ge-ometrical parameters such asvoid volume fractionand voidsize

distribution andthe fracture strength and fracture energy ofthe compositeunderdistinctloadingmodes.

The work is organized into two parts: the first part in

Section 2 refers to an extension of the model developed in (Turteltaub etal.,2018)to explicitlytake voidsintoaccount. This includesan extension ofthe strain and stress power decomposi-tionstoidentifyvoid-relatedterms.Anewversionofacrack-based Hill-Mandelconditionissubsequentlydevelopedandverified.The second partinSection 3 containsa seriesofparametric analyses designedtostudytheinfluenceofdistinctgeometricalfeatures of thevoids.ThemainfindingsaresummarizedinSection4.

2. Microscale formulation with voids 2.1. Preliminaries

A typicalcrosssection ofaplyina fiber-reinforced composite consistsofcarbonorglassfibersembeddedinanepoxymatrixas showninFig.1.In a hierarchicalmultiscaleformulation,a mate-rial point at a macroscale corresponds to the collective behavior ofall microscalematerialpointsinsideavolumeelement suchas theperiodicvolumeelementhighlightedinthefigure.Thevolume element can be subjected to distinct loading conditions (such as nominalmodeI,IIormixedmode)basedonthechosenvaluesof anaverage(macroscopic)strain tensor



¯imposed ontheelement. Forproblems without fracture, periodic boundary conditions are generallypreferredsincethey conform morecloselytotheactual materialbehaviorandallowfasterconvergenceintermsofthesize ofa representative volume element compared to linear displace-mentconditionsoruniformtraction.Forproblemsinvolving frac-ture, periodicboundary conditions can alsobe applied since this allowsmodellingofamacroscopiccrackwithouttheneedto spec-ifyapriorithecrackorientation(see e.g.,(Turteltaub etal., 2018) fordetails).Alternativeboundaryconditionsmaybeused(see,e.g., (Coenen et al., 2012a; 2012b; Boscoet al., 2015; Svenning etal., 2016a; 2016b), but pointwise periodic boundary conditions have

theadvantageofbeingrelativelysimpletoimplementandprovide sufficientflexibilityforthesimulations.

Atypicalmicroscopic volumeelement witha periodiccrackis shown inFig. 2. The volume element alsoshows voids, some of whichmayendupinthepathofacrack(ormayinfactbe a lo-cationwhereacrackoriginates).Thebold lineindicatessegments ofacrack(includingpossibledisconnectedsegments),allofwhich haveanetcontributiontothemacroscopicfractureproperties, par-ticularlytheeffectivefractureenergy.

Inordertodeterminetheeffectivefracturepropertiesofa com-posite, the microscopic volume element is loaded with a given average strain



¯ until complete failure of the specimen. Denote as V =∪nv

k=1V (k) the collection of all voids V (k), k =1,...,n v, in-side thevolume element



anddenoteas



thecollectionof all crackedsegmentsattheendofaquasi-staticprocessthatleadsto complete failure. Using the nomenclature indicated in Fig.2, the periodicboundary-valueproblemthatsimulatesthecracking pro-cessis,intheabsenceofbodyforces,formulatedas

⎧

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎨

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎪

⎩

div

σ

(

x,t

)

=0 xin



\

(



∪V

)

t+

(

x+,t

)

=−t−

₍

_x−_,t

₎

_xon

_

t

(

x,t

)

=tv

₍

_x,t

₎

_xon

_∂

_V u

(

x+l1e1,t

)

− u

(

x,t

)

=l1



¯

(

t

)

e1 t

(

x+l1e1,t

)

=−t

(

x,t

)

xon

∂



3

\

(



∪

∂

V

)

u

(

x+l2e2,t

)

− u

(

x,t

)

=l2



¯

(

t

)

e2 t

(

x+l2e2,t

)

=−t

(

x,t

)

xon

∂



4

\

(



∪

∂

V

)

u±

(

x±+l1e1,t

)

− u±

(

x±,t

)

=l1



¯

(

t

)

e1 t±

(

x±+l1e1,t

)

=−t±

(

x±,t

)

xon

∂



3∩

(



∪

∂

V

)

u±

(

x±+l2e2,t

)

− u±

(

x±,t

)

=l2



¯

(

t

)

e2 t±

(

x±+l2e2,t

)

=−t±

(

x±,t

)

xon

∂



4∩

(



∪

∂

V

)

(1)

where

σ

isthestresstensor,divisthedivergenceoperator,tisthe tractionvectoractingonthecorrespondingsurface(



,

∂



or

∂

V ),

u is the displacement vector and



¯=



¯

(

t

)

corresponds to a pre-scribedmacroscopicstrain tensorapplied onthevolume element thatdrives thedeformationprocess atdifferenttimes t .The trac-tion on the pores is given by a prescribed function tv_{, which is} typically zerounlessthepores are filledwitha fluidunder pres-sure. The set



b_:=

_

\

₍

_

_∪V

₎

_{refers to points in the bulk and} thesuperscripts +and-refer tovalueson oppositessidesofthe surface



asillustratedinthefigure.

Themicroscopicfractureprocessissolvednumericallyusing in-trinsiccohesiveelementsembeddedwithinthebulkelements.This techniqueallowstomodeleach materialphaseseparately aswell astheinterfacesbetweenthefibersandthematrix.Thebehavior of thematrix andthe fibers can be specified witha constitutive modelforthecorrespondingbulkelements,whilethefracture be-haviorismodeledwithtraction-separationrelations(i.e.,cohesive relations see, e.g., Park andPaulino, 2011) which depend on the individual fracture strength andfracture energyofeach phase or interface.

Care must be exercised during the post-processing of results to guarantee that the micro and macro scales quantities are energetically-consistent.Inpracticethismeansthatone hasto ei-therguaranteeaprioriortocheckaposteriorithattheproductof theeffectivemacroscopicquantities(suchasstressandstrainrate) isequaltothevolumeaverageoftheproductofthecorresponding microscopicquantities.Thisprerequisite inhierarchicalmultiscale modelling, known as the scale transition condition, the macro-homogeneity condition or, more commonly, as the Hill-Mandel condition, is an essential requirement that needs to be satis-fiedinanenergetically-consistentmultiscaleformulation(see,e.g.,

Saebetal.,2016 andthereferencestherein).Forproblems involv-ing fracture,the Hill-Mandel conditionrequiresspecial treatment since damagelocalizesin interfaces(Gitman etal.,2007; Nguyen etal.,2010;Unger,2013).Inthepresentwork,a crack-based Hill-Mandelcondition developedin Turteltaub etal.(2018) is further extendedtoincludethepresenceofvoids.Inmanycasesof prac-ticalinterest the voids may be assumedto be traction-free from

S. Turteltaub and G. de Jong / International Journal of Solids and Structures 165 (2019) 63–74 65

Periodic volume elements

One exit/entrance point

Two exit/entrance points

m

f

m

f

One exit/

entrance point

One exit/

entrance point

Two parallel cracks

Three parallel cracks

crack element

crack element

Fig. 1. Periodic volume element in a cross-section of a fiber-reinforced composite. The type of loading is specified through the chosen values ¯of the average (macroscopic) strain tensor. As shown in Turteltaub et al. (2018) , the corresponding localized crack repeats itself periodically but its effective orientation is not geometrically constrained by periodicity and can be interpreted as a single (isolated) crack while parallel periodic cracks may be ignored.

Fig. 2. Nomenclature used in a microstructural volume element . The dashed lines indicate a partition of the domain used to compute integrals in subdomains.

the outset, so they donot directly contribute to the fracture en-ergy(i.e.,tv_{=0 in(1)).However,eveninthetraction-freecase,the} presenceofvoidsneedstobetakenintoaccountwhencomparing quantities such asthe averagestress power per unit volume and the stress powerratefromthe productof theaverage stress and theaveragestrain rate.Tothisend,thestrain isdecomposedinto bulk,crackandvoidcontributionsinordertoseparatetheeffective quantitiesassociatedtothefracture processfromthoserelatedto thebulkprocess.

2.2. Strain decomposition

Theappliedstrain



¯ actingonavolumeelement



isby defi-nitiongivenas ¯



:= 1

|



|

 ∂[un]symds

withthenotation[·]sym indicatingthesymmetricpartofthe ten-sor,nrepresentingtheoutward normalunitvectorasindicatedin

Fig.2, correspondingtothetensorproductand|



|thevolume (areaperunit depth)ofthe volumeelement



.Inorderto iden-tifythepartofthedeformationinthevolumeelement associated tothefractureprocess, considerapartition of



intosubdomains asindicated by dashed lines inFig. 2. The partition is such that the subdomains contain all cracked surfaces



andall void sur-faces

∂

V .Tocoverthewholedomainusingsimply-connected sub-domains, the partition can be complemented by divisions across uncrakedbulkmaterialwherethedisplacementfieldiscontinuous. Asindicated inFig.2,ifthedisplacementfield iscontinuous, the correspondingjumpiszero,(i.e.,[[u]]=0 , with[[

(

·

)

]]=

(

·

)

+₋

₍

_·

₎

− denotingthejumpofaquantity (· )acrossasurfaceandwiththe superscripts+and-referringtothevaluesofthequantityon op-positessidesofthesurface.)Makinguseofthedivergencetheorem ineachsubdomainseparatelyandcollectingallintegralsitfollows thattheappliedstrain



¯canbedecomposedasfollows:

¯



=



b₊

_

f₊

_

v ₍₂₎

where the bulk strain



b_{, the fracture strain}

_

f _{and the “void} strain”



v_{aredefinedas}



b _:= 1

|



|

 b



d

v



= 1₂



∇

u+

(

∇

u

)

T





f _:= 1

|



|

 [[[u]]m]symds



v _:=− 1

|



|

 ∂V [um]symds (3)

withthe superscriptTreferring tothetranspose. Observethat in thedefinition ofthe void strain there is a negativesignthat re-flectsthesignconventionindicated intheinsetinFig.2whereby thenormalvectorm₌mv_{onthevoidsurface}

∂

_{V isdefinedinthe} outwarddirectionfrom thebulk andhenceinwardswithrespect tothevoid.Inthiswayan increase inthevolumeofthevoid,for examplewithuandmorientedinthesamenormaldirectionbut pointinginoppositesenses,resultsina positive normalvoidstrain.

2.3. Stress power decomposition

The stress power P ext_{(power per unit volume)appliedonthe} externalboundary ofthe volume element canbe decomposed in termsofthestresspowerinthebulk, P b_{,therateofworkonthe} cracksurface, P f_{,anda powertermrelatedto thevoids,denoted} as P v_{.Themainingredientrequiredforthestresspower} decompo-sitionisthebalanceoflinearmomentumforaquasi-staticprocess (withoutbodyforce,asindicated in(1)),whichisthenmultiplied bythe velocity field u˙ andintegratedby parts inthe volume el-ement.Similar to the strain decomposition, thestress power de-compositionrequiresrepeateduseofthe divergencetheorem,for whichonecanpartitionthevolumeelement



intosubdomainsas indicatedinFig.2.Asbefore,thedivergencetheoremisused sep-aratelyineachsub-domainandtheintegralsarethenaddedusing thesignconventionsindicatedinFig.2,whichyields

Pext=Pb+Pf− Pv ₍₄₎ where Pext:= 1

|



|

 ∂t· ˙uds Pb _:= 1

|



|

 b

σ

· ˙



d

v

Pf_:= 1

|



|

 t· [[u˙]]ds Pv_:= 1

|



|

 ∂V tv_{· ˙uds. (5)}

Itis in principlepossible to formally view the boundariesof the pores,

∂

V , as an “external” surface with a prescribed “external”

traction but, since it is more natural to view it as an internal surface with a known traction, the term P v _{is accounted for on} therighthandsideofthedecomposition(4).Correspondingly,the term“external” stresspower, P ext_{,is reservedfortherateofwork} done on

∂



.Hence,ifaporeintersects theboundary

∂



,itmay bestill interpretedasan internalboundary.Theterm P b_accounts forthestresspowerintheuncrackedmaterial(bulk)thatoccupies theregion



b_{asindicatedabove.Forthefractureprocess,the} rel-evant termis the rateof work onthe cracked surface, P f_,which does not include the pores. Finally, the term P v_{, which accounts} fortherateofworkdone byanimposedtractiont=tv _actingon thesurfacesofthepores,maybe relevantforsituationswhenthe poresarefilledwith,forexample,apressurizedfluid.Inthatcase aseparate modelisrequiredtodeterminetheappliedtractionon thepores.However,formanystructuralapplicationsofinterest,it may be assumed that this term is negligible compared to other loadsappliedtothestructure.

2.4. Identification of crack, solid and pores contributions

Intheframework ofhierarchicalmultiscaleanalysis,the effec-tive(macroscopic)quantitiesare oftendefinedfromtheoutsetas volume or surface averages of their microscopic counterparts. If the deformation does not localize inside a volume element (e.g., no localized plasticity or cracks appear), then the volume-based Hill-Mandel condition can be satisfied using, for example, peri-odicboundaryconditionsinordertoguaranteeaprioriconsistency across length scales. In that case the product of the averages is equaltotheaverageoftheproducts.

In the context of localized failure, however, the situation is more complexin the sense that the effective quantities in local-izedmodelsmaynotexpresseddirectlyasasurfaceorvolume av-erageoftheirmicroscopiccounterparts.Consequently,inthatcase theHill-Mandel conditionmaynot be satisfieda priori.Asimple solutiontothisissueistoinverttheroleoftheHill-Mandel condi-tion.Usingthisapproach,theeffectivemacroscopicquantitiesare in fact defined(or calibrated) fromthe Hill-Mandel condition. In thatfashion,thetraditionalvolume(orsurface)averagedefinition ofaquantity isreplacedbyan energy-baseddefinitionthat auto-matically guarantees consistencyacross scales. Periodic boundary conditionsmaystillbe usedtosatisfythecorresponding volume-basedHill-Mandelconditionanddefinitionsofeffectivequantities forfieldsassociatedtothebulkmaterialintheneighborhoodofa localizedcrack.

A crack-based Hill-Mandel condition, in a domain containing voids,canbeestablishedfromthedecompositionofthestrainand the stress powerasdeveloped inSections2.2 and2.3.The guid-ing principle in the foregoing analysis is the need to establish a relationbetweenthe rateofwork associatedto the fracture pro-cess,i.e., P f_{asdefinedin(5),andtheproductbetweenaneffective} cohesive traction tf _{andan effective crack opening rate [[u˙]]}f _for ahomogenizedtraction-separationrelationthatmaybeusedata macroscale. Simultaneously,a separaterequirementisto preserve thevolume-baseddefinitionsofeffectivefieldscommonlyusedin multiscale analysisin the absence oflocalization, particularlyfor therateofwork P b_{asindicatedin(5).}

It is convenient to introduce the following notationfor bulk-averagedandcrack-averagedquantities:



(

·

)



:=

_|

_

1

_|

 b

(

·

)

d

v



(

·

)



:=

_|

_

1

_|

 

(

·

)

ds.

Usingperiodic boundaryconditions,itcanbe shownthatthe ex-ternalstresspower P ext _{definedin(5)satisfiesthefollowing} rela-tion:

Pext₌

_

_σ

S. Turteltaub and G. de Jong / International Journal of Solids and Structures 165 (2019) 63–74 67

that corresponds to the classical (volume-based) Hill-Mandel macrohomogeneity condition. Consequently, the effective (macro-scopic)stresstensorcorrespondstothevolumeaverageofits mi-croscopiccounterpart.However, upontheformationofalocalized crack,thetractionvectoronplanesotherthanthecrackitselfmay become discontinuous, hence the stress tensor is discontinuous. The traction across thecrack surfaceduringthe cracking process remains continuous. The purpose of the foregoing analysis is to establisha relationforthetractiononthe crackedsurfaceduring decohesion. Usingthedecompositions(2)(intime-rateform) and

(4)intherelationaboveyields

Pb_+Pf_{− P}v₌

_

_σ

·





˙ +



σ

· ˙



f+



σ

· ˙



v.

Consequently,fromthestrainandpowerdecompositions,itis nat-uraltoidentifythreeseparate relationsassociatedtoscale transi-tionsasfollows:

Rateofworkin/on Microscaleterm Proposedmacroscaleequivalent

Bulk Pb

_

σ

·





˙ 

Cracksurfaces Pf

_

σ

· ˙



f

Void/poresurfaces Pv ₋

_

_σ

· ˙



v

(6)

Observe that thenegative signin theexpression for thevoids is a consequenceofthesignconventionwhereby thenormalvector pointsoutwardsfromthesolidandintothevoid.Itshouldbe em-phasized that theaforementioned identifications are assumptions thatneedtobeverifiedaposteriori.

2.5. A crack-based Hill-Mandel condition

From amodellingpointofview,poresmaybehandledintwo possibleways:(i)treattheporesasentitiesseparatefromthebulk and the crack or (ii) combine pores with the effective behavior of the bulk and/or the crack.For example, pores that are in the path ofacrackmaybe combinedwiththecrackitselfwhile iso-latedporesmaybe combinedwiththe bulkmaterial.However, if the traction prescribedon the poresurfacesis not zero, combin-ing the pores with the crack may cause modelling problems in the formof a residual cohesive traction that is connectedto the loadinganddoesnotrepresenttheactualfracturebehaviorofthe material.Toresolvethisissue,amorerobust approachisto com-bine all thepores withthebulk behaviorregardlessoftheir loca-tion(i.e.,whethertheyareonthecrackpathornot).Poresinduce stress concentrationsandfurther facilitatecracking,which gener-allyaffects thefracture strength andenergy.Hence, althoughthe poresare formallynotconsidered aspartofthecrack,their pres-encehasastronginfluenceontheeffectivetraction-separation re-lations.In linewiththisconsideration,unlessexplicitlyindicated, the “bulk” refers henceforth to the combined solid and porous parts while the “crack” refers to the surfaces created duringthe cracking process but excluding pre-exisiting pores. Consequently, the crack-based Hill-Mandel condition is treated separately from thecombinedsolidandporousparts.

AsdiscussedinNguyenetal.(2010)andTurteltaubetal.(2018), to prevent RVE convergence problems related to non-constant surface-to-volumeratiosforvolumeelementsofdistinctsize,itis important to averagefracture properties on the cracksurface in-steadthanonthevolumeelement.Tothisend,multiplythe crack-basedstress power P f _{in(3)andthe correspondingterm}

_

_σ

· ˙



f

(as shownin(6)) by thefactor|



|/|



| to get,inview of(5), the followingcorrespondencerelation:



t· [[u˙]]



_↔



σ

·

[[[u˙]]m]sym

. (7)

The symbol ↔ isused in(7) toemphasize that equalityisnot a priori guaranteedbut, rather, itis aproposed scale transition re-quirement.

2.6. Effective fracture behavior

For a cohesive-zone modelling approach at the macroscopic level,the model should involve two ingredients: (i) a nucleation criterion(i.e.,underwhatloadingconditionsandinwhich orienta-tionacrackisfirst detectedatthemacroscopiclengthscale)and (ii)anevolutionmodel(i.e.,howdoesthecohesivetractionevolves withcrackopening).Althoughaninelasticprocessshouldnormally be expressed in rate form to account for path-dependency, it is customaryinthecohesivezone frameworktomodelthecracking process using directly a relation between an (effective) cohesive tractiontf_{andaneffectivecrackopening[[u]]}f_{.Inaddition,dueto} thescaletransition,themacroscopicdescriptionrequiresan effec-tivecracklength|



f_{|associatedwithamacroscopiccrackwithina} volumeelementasillustratedinFig.3.

Animportantissue isthedefinition oftheeffective quantities |



f_{|, t}f _and[[u]]f_{.Observethat,incombinationwith(3),theright} handsideofrelation(7),whichshouldrefertotheeffective prop-erties,involvestensorialquantities(i.e.,



σ

_ and



f_)andnot vec-tors(i.e.,tf _and[[u]]f_{) ascommonly-usedina traction-separation} relation.Therelationbetweenthesequantitiesisnottrivialsinceit involvesitselfascaletransition(e.g.,



t



₌



σ

n



isgenerally differ-entthan



σ

n



andasimilar observationapplies forthe fracture straintensorandcrackopeningvector).

From a geometrical point of view it is natural to define the effective crack opening rate as the surface average of its micro-scopiccounterpart,scaledwiththeratioofeffectivetoactualcrack length,i.e.,

[[u]]f:=



|



|



f



[[u]]



 (8)

where the effective crack length |



f_{| is obtained from the} ge-ometry of a volume element using the approach proposed in

Turteltaubetal.(2018),i.e.,



_

f



_:=

⎧

⎨

⎩



_

f min



ifr≥ rmax



_

f max



ifr<rmax (9)

Fig. 3. Illustration of a sub-ply volume element with a microscopic crack (right) and equivalent macroscale crack (left) with effective normal vector m f , effective co-

hesive traction t f and effective crack opening [[ u ]] f . Observe that here the notion

of macroscale refers to a larger scale compared to the microscale but still within a single ply in a laminate.

where,usingthenotationshowninFig.2,



_

f min



:=min



l1



n2· mf



, l2



n1· mf





_

f max



:=max



l1



n2· mf



, l2



n1· mf



r:=





f max





_

f min



(10)

andwiththeeffectivecracknormalbeingdefinedas

mf_:=

_

_m

_

. (11)

Theprevious expression isbasedon thenotion thatthe effective cracklength|



f_{|shouldcorrespondtothelengthofastraight} seg-mentcrossingperiodicallythevolumeelement inadirection per-pendiculartomf_{.The value r representsthe numberofcrossings} ofastraightperiodiccrackalongaperiodic l 1× l2 rectangular do-main(eitherbottomtotoporleft toright).Ingeneral,however,it isnotpossibletoconstructastraightperiodiccrackwiththegiven vectormf _{(i.e., r isnotanintegerunlesstheorientationofm}f sat-isfiesa specificgeometricalrelationwith l 1 and l 2). Consequently, |



f_{| isinterpreted asan approximation.Further, thecut-off value}

r max is introduced to handle near vertical (or near horizontal) pe-riodiccracks.Inpracticeamaximumvalueof r max=5hasproven tobesufficienttodistinguishbetweenvertical(orhorizontal)and nearvertical (ornearhorizontal)cracks.It isobservedthat,since themicroscalecrackmaycontainmultiplebranches(bifurcations) and/orisolatedsegments,caremustbeexercisedtopreserve con-sistency in the sign conventions of



+ and



− when perform-ingthe integration on



.Further, thisalso impliesthat the ratio |



|/|



f_{| may deviate fromunity, withvalues smaller than 1 also} beingpossibleforcracksthatcontainvoidsintheirpaths.

As indicatedabove,itisconvenienttotie thedefinitionofthe effectivecohesivetractiontotherelation(7),i.e.,tf_canbe implic-itlydefinedsuchthat



_

f



_tf_{· [[u˙]]}f_:=

_|

_

_|

_

_σ

·

[[[u˙]]m]sym

. (12)

Correspondingly, in view of (7), (8) and (12), tf _{is defined such} thatthemacroscopicdescriptionoftherateofworkcoincideswith the microscopic one. One limitation withthe previous definition isthat thepart ofthe tractionthat is perpendicular tothe crack openingrateremains undetermined(i.e.,the formulawouldonly specifytheprojection oftf _{into [[u]]}f_{). Furthermore,froman} im-plementationpointofview,thecalculationoftf_{from(12)maybe} pronetoerrors dueto inaccuraciesin thecomputationofthe in-stantaneousvalues ofthe local crack openingrates. As indicated inTurteltaubetal.(2018),an alternativeapproachistoproposea generalformfor tf _{witha calibration parameter that isthen} ad-justed to satisfy relation (12). Adopting that approach, the pro-posedapproximation for theeffective cohesive traction is as fol-lows:

tf_≈

_α

_tf

+

(

1−

α)

tf (13)

where

α

is a calibration parameterand thecrack-based cohesive tractiontf

 andvolume-basedcohesive tractiontf aredefinedas

tf

:=



t



 tf :=



σ

mf. (14)

The (phenomenological) parameter

α

is taken here as a time-independent quantity that is meant to approximately enforce

(12)ina time-averaged fashion (i.e.,throughouttheentire crack-ingprocess),thusitshould notbeusedtoenforcethescale tran-sitionrelationateachtimeseparately,whichwouldrequirea vari-ableparameter.Toemphasizetherationaleofthisapproach,inthe presentwork(12)isviewedasadefinitionoftheeffectivetraction whereas(13)istakenasanapproximationusedtoenforcethe def-inition. Inview of (8),(11), (12),(13), and(14),the parameter

α

canbeobtainedthroughthepost-processingofsimulationdataof afractureprocessfrom t =0(initialstate)to t = T (fully-damaged state),fromthefollowingrelation:

α

= T 0

(



t



−



σ





m





)

·



[[u˙]]



dt T 0



σ

·

(



[[u˙]]m



−



[[u˙]]







m





)

dt .

After establishing consistency in energy between micro and macroscales, the next step is to determine the actual effective properties(i.e.,convergedquantitieswithrespect tosamplesize). To this end, microscopic volume elementsof increasing size are subjectedtothesamenominalloadtodetermine(numerically)the minimumsizeoftherepresentativevolumeelement (RVE),which corresponds tothesmallestone forwhich thefractureproperties areno longerdependent uponthesize ofthevolumeelement to within aprescribed tolerance(see e.g. Kanitetal.(2003)). Inthe context of localization, one may refer to the equivalentcrack as a representative surface element (RSE) embeddedinan RVE.This analysisrequiresmultiplesamplesofvolumeelementsofthesame sizeto compute themeanofall realizations,althoughinpractice itmayberequiredtousearelatively largeconvergencetolerance and/or to limit the number of samples due to the overall com-putational cost. This is because in each convergence analysisfor thevolumeelementsize,thenumericalresultsmustalsoconverge, thusaseparatemeshrefinementconvergencetestmustbecarried outforeachvolumeelement.Thisconvergencetestisparticularly important in the context of fracture since there is a process of localization(i.e., theformation ofa crack) whichrequires a non-traditionalconvergenceanalysis(seee.g.Turteltaubetal.,2018for details).Once the minimumacceptable size ofthe RVE hasbeen established,parametricanalysescanbeconductedtostudythe in-fluence of physicaland geometrical parameters associated to de-fectsasexplainedinthenextsection.

3. Simulations with volume elements containing voids 3.1. Overview of simulation cases

Inordertostudytheeffectofsub-plyporosities,severalcases areanalyzedusingdistincttypesofsamplesasillustratedinFig.4. Forsimplicity,aplanestrain frameworkisadoptedassuming that inthe out-of-plane direction(longitudinal directionofthefibers) thedimensionsofthesamplearelargecomparedtothesizeofthe computational domain(cross-sectionin theinteriorofa sample). Porositiesof0%,1%,2%,4%,6%and8%wereconsidered(measured perunitcross-sectionalarea).Foreachporosity,twodistincttypes of voids are analyzed, namely matrix voids and interfiber voids. In addition, forthe matrix voids, two typesof cases are consid-eredseparately,namelyone whereall voids havethesame char-acteristic size and one with variable size and aspect ratio using randomly-orientedellipses.Thesetupconsistingofmatrixvoidsof constant size isseen asa baselinefor agiven porosity (withthe caseofno porosity beingtheoverall benchmark tostudythe ef-fectofporosity).Caseswithinterfibervoids(withinterfibervoids up making up about 1/4 of the voids and the restbeing matrix voids),areusedtogagetheeffectofthetypeofmicroscopicvoid. Othereffectsarestudiedbyvaryingthesizeoftheindividualpores andtheiraspectratiowhilekeepingtheporosityfixed,asindicated above.

Since theeffectofdefects isexpectedto depend onthe load-ing conditions, two representative loading caseswere considered foreach configurationasindicatedin Table1.The twocases, ex-pressedintermsoftheapplied macroscopicstraintensor



¯, with

γ

>0 acting as a strain parameter, induce cracking in a cross-sectionperpendiculartothefiberdirection(i.e.,lateralcracking).

Toaccount forstatisticalvariationsinthemicrostructure,eight randomrealizationsareconsideredforeachporosity,foreachvoid

S. Turteltaub and G. de Jong / International Journal of Solids and Structures 165 (2019) 63–74 69

Fig. 4. Illustration of parametric study for defects: for porosities from 0 to 8% distinct configurations are analyzed in terms of void type, shape and size.

Table 1

Load cases expressed in terms of the applied macroscopic in-plane strain ten- sor ¯, with γ> 0 a strain parameter. Observe that the pure shear is applied

at an angle of 45 ° with respect to the coordinate system shown in Fig. 2 .

Load cases Applied in-plane strain

1. Laterally-constrained uniaxial extension

 ¯ 11 ¯12 ¯ 21 ¯22  =  γ 0 0 0  2. Pure Shear  ¯ 11 ¯12 ¯ 21 ¯22  =  γ 0 0 −γ 

typeandshapeandforeachloadcase.Toachievethetarget poros-ityinasample,adiscretenumberofporeswereusedrangingfrom typically6poresfor1%porositytoabout50poresat8% porosity, withanaveragecharacteristicsizeofabout1.7μm.Ingeneralmore realizationsarerequiredtoestablishaccuratestatisticaldata, how-everthenumberislimitedduetothe overallcomputationalcost. To further control the overall computational cost andin view of the largenumberofcasesconsidered above,the material proper-tieswerekeptfixedinall simulationsandallsampleshadafixed fiber volume fraction of(approximately) 50%, withthe rest (ma-trixandvoids)addinguptotheremaining50%.Thematerial prop-ertiesshowninTable 2arerepresentativeofcomposite materials

typicallyencounteredintheaerospaceindustry,anditisexpected thatthefindingsregardingtheeffectofmicrovoidswouldbe sim-ilarforothercompositeswithcomparableproperties.

Due to the large number of simulations, a dedicated python script was developed to generate the samples, generate the meshes,createtheinput filesforthefiniteelement analyses,run thesimulationsinparallel,executethepost-processingand gener-atetheresults. Thesimulations were carriedout usingthe Finite ElementpackageAbaqus(2018)usinglinearplanestrainelements for the bulk (matrix and fiber) and cohesive elementswith lin-earsofteningtocapturethefractureprocess.Unstructuredmeshes were produced using the mesh generator Gmsh (2018) with the goalto minimizethe meshdependency ofthecrack pattern.The cohesiveelementswere embeddedbetweenalledges ofthebulk elements (including fibers). Mesh refinement analysis was sys-tematically carried out for all samples of all sizes until a con-vergedmeshwasidentified withina giventolerance.The conver-genceanalysistoidentifythesizeoftherepresentativevolume el-ement(or,morespecifically,therepresentativesurfaceelementfor an equivalentcrack)wasdone withthe corresponding converged mesh foreach sample size. It wasestablished that fora volume element of75 μm _{× 75} μm themean response waswithin the toleranceofa(discrete)deviationoftheprevioussizevolume

ele-Table 2

Material properties of a typical composite system used in aerospace applications. Elastic properties are taken from HexTow IM7 fibers and CYCOM 5230-2 matrix. Fracture properties are not measured but are assumed to be representative (observe that the fiber-matrix interface is the weakest). Unless otherwise indicated and where applicable, the elastic and fracture properties for the matrix and the fiber are associated to (transversely-)isotropic models.

Property Units Fiber Matrix Fiber-matrix interface Modulus of elasticity GPa 19 (Transverse) 3.52 -

Poisson’s ratio - 0.23 (Transverse) 0.35 -

Fracture strength MPa 100 50 25

Laterally-constrained extension (Mode I) Pure shear (Mode II)

Fig. 5. Typical simulations for the case of 8% porosity under loading case 1 (laterally-constrained uniaxial extension shown on the left) and loading case 2 (pure shear shown on the right). Blue segments indicate crack opening in the matrix, red in the fiber-matrix interface and the rest in the matrix voids and interfiber voids. The effective cracks for loading cases 1 and 2 appeared (on average) as opening modes I and II, i.e., as normal and tangential openings, respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

mentsize62.5μm _{× 62.5}μmwithanominalmeshsize param-eterof1μm,althoughalocallyrefinedmeshwasusedtocapture thedetailsofthevoids.Typicalsimulationsforloadingcases1and 2areshowninFig.5foronerealizationof8%porositywithmatrix voidsandinterfibervoids.

Inmostcasesthecrackpathcontainsportionsthroughthe ma-trix(shown as blue segments in color version), through the rel-ativelyweak matrix-fiber interface (shown asred segments) and throughmatrixandinterfibervoids.

Due toitshigherfracture strengthandenergy,no crackswere observedinthefibers.Further,sincethefibersdidnotcompletely restrict the path of the cracks, an effective crack opening mode typeI (normal)wasobserved inthesimulations forloading case 1while an effectivecrack openingmode type II (tangential) was triggered in the simulations for loading case 2. Correspondingly, theloadingcases1andmaybenominallyinterpretedasfracture modesIandII, respectively,whenreferredtoa crosssection per-pendiculartothefiberdirection.Ingeneral,themodeIsimulations typicallyrequiredtheleastcomputationaleffort,whereas modeII requiredsmalltime stepsto convergeandhencealarger compu-tationaltime. Forconciseness,resultsofthe convergenceanalysis areomittedherebutsomedetailscanbefoundindeJong(2018).

3.2. Matrix voids of constant size

The first set of simulations is carried out for matrix voids of constant size, which will be used subsequently as a benchmark

fora fixedporosity.The effectivetraction-separation relationsfor porositiesrangingfrom0%to8% areshowninFig.6a under con-straineduniaxialextension (nominalmodeI)intermsofthe nor-malcomponent t f

n ofthecohesive tractionvectortf asa function

ofthenormalcomponent[[u n]]foftheeffectivecrackopening

vec-tor[[u]]f_{.Similarly,theeffectivetraction-separationrelationsunder} pureshear(nominalmodeII)areshowninFig.6bintermsofthe tangentialcohesivetractioncomponent t f

sasafunctionofthe

tan-gentialeffectivecrackopeningcomponent[[u s]]f.Thelinesplotted

inthe figures representtheaverageof eightrealizations foreach porosity,while the shaded regions in the graphics correspond to the(discrete)standarddeviation.Forclaritythestandarddeviation isonlyshownonlyfortheporosityof8%.Forloadingcase1 (nom-inalmodeI),the tangentialcomponentwasrelativelysmall com-paredtothenormaloneand,conversely,forloadingcase2 (nomi-nalmodeII),thenormalcomponentwasrelativelysmallcompared tothetangentialcomponent.

Microscalefailure,measuredbymonitoringindividualcohesive elements,initiatedinplacesintrinsicallyweakandpronetostress concentrationsforbothloadingcases.However,microscaledamage doesnotimmediatelytranslateintomacroscopicfailure.Indeed,as itis commoninmultiscaleanalysisoffracture, theonsetof frac-tureat the microscaletypically occurs duringthe early stagesof theloadingbutcan onlybe detectedatthemacroscaleonce suf-ficient damage has accumulated(see, e.g., Suikerand Turteltaub, 2007a;SuikerandTurteltaub,2007b;Farleetal.,2018;Hilleetal., 2011 for examples in distinct material systems). In the current

Fig. 6. Effective traction-separation relation for porosities varying from 0 to 8% for the case of matrix pores of constant size for (a) constrained uniaxial extension (nominal mode I) and (b) pure shear (nominal mode II). See text for further details.

S. Turteltaub and G. de Jong / International Journal of Solids and Structures 165 (2019) 63–74 71

Fig. 7. Effective traction-separation relation for porosities varying from 0 to 8% for the case of matrix pores of variable size for (a) constrained uniaxial extension (mode I) and (b) pure shear (nominal mode II). See text for further details.

context, the first detectable onset of fracture at the macroscale is, by definition, measured by the peak value in the traction-separation relation, which isthe (transverse) macroscale fracture strengthoftheplyunderthecorrespondingloadingcase.

Oncetherepresentativevolumeelementiscompletelydamaged it can no longer transmit loads,in which casethe energy dissi-pated due to fracture in the representative surface element cor-respondstothemacroscopicfracture energy.Asmaybeobserved from Fig. 6a, the fracture strength for mode Iremains relatively constant fordistinctporositieswhile thefractureenergy(visually correspondingtotheareaunderthetraction-separationcurve), de-creases for increasing porosity. For mode II, as may be seen in

Fig.6b,boththefracturestrengthandthefractureenergydecrease asa functionof porosity.Observethatin Fig.6bsome curvesdo not terminate atzero traction; this is dueto the fact that some oftherealizationshadconvergenceissuesinshearwhenthecrack wasalmostfully-formed.ThecurvesinFig.6b,whichareobtained asan averageofmultiplerealizations,onlyreflecttheportion for whichallrealizationsconverged.Forthepurposeofcomputingthe effectivefractureenergy,whichrequiresaresponseuntilcomplete damage, an interpolation to zero traction wasused for the non-converged curvesandaveraged withtheconverged curvesofthe sameporosity.The effectivepropertiesobtainedthrough postpro-cessing of the simulations, particularly the fracture strength, are also strongly influenced by the fiber-matrix interface properties since thisisthe weakest componentinthe system(see Table2). Aquantitativeanalysisofthedependencyonporosityispresented belowincombinationwiththeresultsfromtwoothersetsof sim-ulations,namelymatrixvoidsofvariablesizeandvoidsofdistinct types.

3.3. Matrix voids of variable size and aspect ratio

Thesecondsetofsimulationsiscarriedoutformatrixvoidsof variable size andaspect ratio.The aspect ratiosvaried from 1to 1.4 andthe characterisitc sizes were derived froma normal dis-tribution around the average characterisitic size. Due to the lim-ited numbers of pores, however, only a few representative sizes were considered in each simulations ranging from 0.17 to 1.83 times the average characteristic size. As in the previous bench-markcaseofconstantsize,resultsareshowninthissectionforthe twoloadingcasesandvariousporositiesinFig.7,rangingfrom0% to8%. Inparticular,Fig.7acorrespondsto thenormalcomponent

t f

n of the cohesive traction vector tf asa function ofthe normal

component[[u n]]f ofthe effectivecrackopening vector[[u]]f

un-derconstraineduniaxialextension(nominalmodeI),whileFig.7b showstheeffectivetangentialcohesive tractioncomponent t f

sasa

functionofthetangentialeffectivecrackopeningcomponent[[u s]]f

under pure shear (nominal mode II). As before, the lines indi-catetheaverageresponse,overeightrandomrealizationsforeach porosity,while theshaded areas illustratethe standard deviation (shownonlyfortheporosityof8%).

Similarqualitativeobservationscanbedrawninthiscase com-paredto the previous set ofsimulations, namelya weak depen-dencyofthefracturestrengthon porosityunderconstrained uni-axial extension, a decrease of the effective strength under pure shear withincreasing porosity and, forboth loading cases,a de-creaseinfractureenergywithincreasingporosity.

3.4. Matrix and interfiber voids

The third setof simulations is carriedout for interfiber voids together with matrix voids of fixed size. As before, results are shownin Fig. 8for the two loading casesand various porosities (0%to8%).Forconstraineduniaxialextensionthenormal compo-nent t f

n ofthecohesive traction isshowninFig. 8a asa function

of the normal component of the effective crack opening, [[u n]]f,

whereasforpureshearthetangentialcohesivetractioncomponent

t f

sisshownasafunctionofthetangentialeffectivecrackopening

component[[u s]]f in Fig.8b. Similar to theprevious sets of

sim-ulations, the solid lines are obtained from the average of eight realizations while the shaded regions indicate the discrete stan-dard deviation. In this caseit is worth noting that, dueto com-putational andgeometrical limitations, only a limited amount of voidscorrespondtointerfibervoidssincetheseregionsareformed bygroupsofthree orfourfibersthat areincloseproximity (see, e.g., Fig. 9). Although the net average distance between fibers is not greatly affected(i.e., negligiblefiber clustering),these factors need tobe takeninto account in orderto gage the neteffect of interfibervoids.Acomparativeanalysisofthesimulationresultsof

Figs.6–8ispresentedinthenextsection.

3.5. Comparative analysis

Inordertoquantifytheinfluenceofthegeometricalfeaturesof themicrostructureontheeffectivefractureproperties,itis conve-nienttosummarizealltheresultsforeachloadingcaseseparately.

Fig. 8. Effective traction-separation relation for porosities varying from 0 to 8% for the case of interfiber voids and matrix pores of fixed size for (a) constrained uniaxial extension (mode I) and (b) pure shear (nominal mode II). See text for further details.

Fig. 9. Normal component of the effective fracture strength under uniaxial exten- sion (mode I) as a function of porosity for different types and distributions of de- fects. The tangential component of the effective strength is negligible.

3.5.1. Comparative analysis for constrained uniaxial extension (mode I)

The normalcomponent t nf,croftheeffective(macroscopic)

frac-turestrengthunderuniaxialextensionisshowninFig.9asa func-tionofporosity forthethreesetsofvoidsanalyzed(matrixvoids of constant size, matrix voids of variable size and a mixture of interfiber voids and matrix voids). In all simulations, the cracks in the numerical samples open in a nominal mode I. The error bars ateach porosity analyzed indicate thediscrete standard de-viationfromthemean.Tovisualizethedependencyofthefracture strength on porosity, the shaded regions indicate the upper and lowervalues ofthe standard deviationwhile themeans foreach set are connected with dashed lines. In addition, a linear fit for themeansisshownasasolidline.

From Fig. 9, it may be observed that the effective fracture strength isdominated by the weakest componentin the system, namelythefiber-matrixinterfaces,withaneffectivevalue slightly above25MPa(see Table 2).Thisvalue, interestingly,remains rel-ativelyconstantas afunction ofporosity.Although theinterfiber voids casehas a slightlylower effective strength, the differences arenotsignificantcomparedtothematrixvoids(witheither con-stantorvariablevoidsize).Itisrelevanttorecallthattheamount ofvoidsdirectlyconnected tointerfibervoids isabout1/4of the porosity,whiletherestiscomposedofmatrixvoids,hencethe ef-fectduetointerfibervoidsiscoupledtotheeffectofmatrixvoids.

Fig. 10. Effective fracture energy G f

I under uniaxial extension (mode I) as a function

of porosity for different types and distributions of defects.

A directcomparison ofpure interfibers voids versus pure matrix voids was not carried out due to the limited range of porosity achievable withpure interfiber voids(up to2%). Forthe0% case (novoids),allthreecasesshouldinprinciplestartfromacommon point.However,thedeviationsshownintheresultsatzero poros-itybetweenthethree typesofvoids analyzedcan beascribed to thefactthatthesewereobtainedfromdistinctsetsofrealizations, hencethedistinctvaluesreflecttheinfluenceofthefiber distribu-tionsandnotofthevoids.Ingeneral, thesimulationresults indi-catethat, contrarytoexpectations,the fracture strengthin mode Iisrelativelyinsensitivetodefectsuptotherangeofvalues ana-lyzed.

The effective fracture energy G f

I for the nominal mode I is

shown in Fig. 10 as a function of porosity for the three sets of voids analyzed.From thisfigure,it can be seen that the fracture energydecreases significantly withincreasing porosity. Thevalue forzeroporosityisstronglyinfluencedbytheweakestcomponent, namelythe fiber-matrixinterfaces (see Table2). The fracture en-ergydecreasesbyabout50%(basedonmeanvalues)asthe poros-ityincreasesto 8%,hence ithasa strongdependency ondefects. The mean values ofthe fracture energy for the three cases con-sideredaresimilarandsoistherateofdecreaseasshowninthe solidlines(linearfit).Theseresultsindicatethatthemain param-eteraffectingthedecreaseinfractureenergyistheporosity,while thesizeandtypeofthevoidsonlyplayasecondaryrole.

S. Turteltaub and G. de Jong / International Journal of Solids and Structures 165 (2019) 63–74 73

Fig. 11. Tangential component of the effective fracture strength under shear defor- mation (mode II) as a function of porosity for different types and distributions of defects. The normal component of the effective strength is negligible.

3.5.2. Comparative analysis for pure shear loading (mode II)

Forapuresheardeformation,thetangentialcomponent t sf,crof

theeffective(macroscopic)fracturestrengthisshowninFig.11as afunctionofporosityforthethreesetsofvoidsanalyzed.The er-ror bars ateach porosity analyzed indicate the discrete standard deviation from the mean of the fracture strength. In all simula-tions,thecracksinthenumericalsamplesopeninanominalmode II.SimilartothemodeIcase,thefracturestrengthforthemodeII caseatzeroporosity(novoids)isinfluencedbytheweakestphase (fiber-matrix interface). The strength is about 20% to 30% higher thaninmodeIforthedistinctrealizations.Thesimulationsatzero porositywere carriedout fromdistinct setsofmicrostructures for eachvoidconfiguration,hencethedistinctvaluesreflectthe influ-enceofthefiberdistribution.From theresultsinFig.11itcanbe observed that the mode II strength decreases approximately lin-early by about 15%to 30% asthe porosity varies from0% to 8%. Therateofdecreaseappearsslightlylargerfortheinterfibercase, althoughthedifferenceismarginallysignificant.

The effective fracture energy G c

II under mode II is shown in

Fig. 12 as a function of the porosity forthe different types and distributionsofvoids.SimilartothemodeIcase,themodeII frac-tureenergyforthezeroporositycaseisdependentuponthe fiber-matrix interface properties (weakest phase). The mean value for zero porosity is about50% higherthan the zeroporosity mode I fracture energy. The fracture energyin modeII decreases signifi-cantlyby about60% astheporosity increasesto 8%,a somewhat larger drop compared to the mode I behavior. Further, the data suggeststhatthedecreaseisbetterapproximatedwithaquadratic function that starts to saturate at the end (i.e., at 8% porosity).

Fig. 12. Effective fracture energy G f

II under shear deformation (mode II) as a func-

tion of porosity for different types and distributions of defects.

Observethattheinterfibercaseappearstohaveananomalyat1% withaslightly highermeanvalue compared tothe novoidcase, but this difference is not sufficiently significant. It may also be interpreted as a relatively constant mode II fracture strength for smallvaluesoftheporosity,buttheoveralltrendisthatthe frac-tureenergydecreaseswithporosityinall cases.As inthecaseof modeI,thedifferencesbetweenthetypesandsizesofvoidsonly appeartoplayasecondaryrolewhilethemainsignificant param-eteristheporosity.

Experimentalresultsfortransversefracturestrengthasa func-tionof porosity withthe sameset-up asthe onesanalyzed here donot appear to be readily available in theliterature. Neverthe-less,itisworth pointingout thatexperimental measurementsfor out-of-planetensilestrengthofalaminate aswellasinterlaminar shear strength are qualitatively similar, with reductionsof about 50% whenthe porosity increasesfrom zeroto 7% asreportedby

Bowles andFrimpong (1992); Gurdal etal. (1991); Jeong(1997);

Costaetal.(2001).Althoughcomparisonsbetweenthe aforemen-tionedexperimentsandthepresentresultsshouldnotbe donein aone-to-onefashion, theexperimental dataprovides indirect ev-idencethat the orderofmagnitude ofthepresentsimulation re-sultsisrepresentativefortheactualmaterialsystemanalyzed.

4. Concluding remarks

Anewmultiscaleaveragingmethodhasbeendevelopedto ex-plicitlystudytheinfluenceofmicrovoidsintheply-levelbehavior ofa unidirectional laminate. The methodcan be readilyused for moregeneralcompositematerials orformicromechanicalanalysis offractureinothermaterials(metals,concrete,etc.).Thekey ingre-dientinthemethodistheseparationofthecontributionofvoids, cracksandbulkandtheenforcementoftheHill-mandelcondition fora representativesurfaceelement(effectivecrack).The method wasappliedtostudytheinfluenceofvoidsonthetransverse frac-turestrengthandfractureenergyoffiber-reinfocedunidirectional composites.Themainconclusionsfromtheanalysisofmicrovoids inafiber-reinforcedcompositeareasfollows:

• The effectivemode Ifracturestrengthremains nearlyconstant withincreasing porosity whilethere isa largelinear decrease intheeffectivemodeIfractureenergy.

• There is a significant decrease in both the effective mode II fracture strength (linearly) and energy (quadratically) for in-creasingporosity

• The porosity is the main parameter influencing the effective fractureproperties.

• Microscopic voidsize andtype only play a secondary role in affectingthefractureproperties.

The findings in the present study are limited to transverse propertiesduetocomputationalissues,butthemethodologymay be extended to coupled three-dimensional cases. Nevertheless, it isexpectedthatsimulationsinvolvingfiberandmatrixcracksmay posecomputationalchallengesduetothelargedifferencesin frac-ture properties. Further, other loading cases, in particular com-pressiveloading,mayrequirespecial treatmentduetothe unsta-ble fracture mechanisms involved such askink bands. Neverthe-less,thesearetechnologicallyrelevantcasesthatrequireattention particularlysince defectsplay a criticalrole incompressive load-ing(Kosmann etal.,2015), andthepresentmethodologymaybe usedtostudythesecasesaswell.

References

Abaqus, 2018. Finite Element Analysis. www.3ds.com/products-services/simulia/ products/abaqus/ . Accessed: 2019-1-21.

Arteiro, A. , Catalanotti, G. , Melro, A.R. , Linde, P. , Camanho, P.P. , 2015. Micro-mechani- cal analysis of the effect of ply thickness on the transverse compressive strength of polymer composites. Compos. Part A 79, 127–137 .

Bosco, E. , Kouznetsova, V.G. , Geers, M.G.D. , 2015. Multi-scale computational homog- enization-localization for propagating discontinuities using X-FEM. Int. J. Numer. Methods Eng. 102 (3–4, SI), 496–527 .

Bowles, K. , Frimpong, S. , 1992. Void effects on the interlaminar shear-strength of unidirectional graphite-fiber-reinforced composites. J. Compos. Mater. 26 (10), 1487–1509 .

Canal, L.P. , Segurado, J. , Llorca, J. , 2009. Failure surface of epoxy-modified fiber-rein- forced composites under transverse tension and out-of-plane shear. Int. J. Solids Struct. 46 (11–12), 2265–2274 .

Coenen, E.W.C. , Kouznetsova, V.G. , Bosco, E. , Geers, M.G.D. , 2012a. A multi-scale approach to bridge microscale damage and macroscale failure: a nested com- putational homogenization-localization framework. Int. J. Fract. 178 (1–2, SI), 157–178 .

Coenen, E.W.C. , Kouznetsova, V.G. , Geers, M.G.D. , 2012b. Novel boundary conditions for strain localization analyses in microstructural volume elements. Int. J. Nu- mer. Methods Eng. 90 (1), 1–21 .

Costa, M.L. , de Almeida, S.F.M. , Rezende, M.C. , 2001. The influence of porosity on the interlaminar shear strength of carbon/epoxy and carbon/bismaleimide fabric laminates. Compos. Sci. Technol. 61 (14), 2101–2108 .

Farle, A.-S. , Krishnasamy, J. , Turteltaub, S. , Kwakernaak, C. , van der Zwaag, S. , Sloof, W.G. , 2018. Determination of fracture strength and fracture energy of (metallo-) ceramics by a wedge loading methodology and corresponding cohe- sive zone-based finite element analysis. Eng. Fract. Mech. 196, 56–70 . Gitman, I. , Askes, H. , Sluys, L. , 2007. Representative volume: existence and size de-

termination. Eng. Fract. Mech. 74 (16), 2518–2534 .

Gmsh, 2018. Finite Element Mesh Generator. www.gmsh.info . Accessed: 2019-1-21.

Gurdal, Z. , Tomasino, A. , Biggers, S. , 1991. Effects of processing induced defects on laminate response - interlaminar tensile-strength. SAMPE J. 27 (4), 39–49 . Hamidi, Y.K. , Aktas, L. , Altan, M.C. , 2004. Formation of microscopic voids in resin

transfer molded composites. J. Eng. Mater. Technol. 126 (4), 420–426 . Hettich, T. , Hund, A. , Ramm, E. , 2008. Modeling of failure in composites by X-FEM

and level sets within a multiscale framework. Comput. Methods Appl. Mech. Eng. Enriched Simulation Methods and Related Topics 197 (5), 414–424 . Hille, T.S. , Turteltaub, S. , Suiker, A.S.J. , 2011. Oxide growth and damage evolution in

thermal barrier coatings. Eng. Fract. Mech. 78 (10), 2139–2152 .

Jeong, H. , 1997. Effects of voids on the mechanical strength and ultrasonic attenua- tion of laminated composites. J. Compos. Mater. 31 (3), 276–292 .

de Jong, G. , 2018. Multiscale Modeling of the Effect of Sub-Ply Voids on the Failure of Composites: A progressive Failure Model. Delft University of Technology, Delft Msc thesis .

Kanit, T. , Forest, S. , Galliet, I. , Mounoury, V. , Jeulin, D. , 2003. Determination of the size of the representative volume element for random composites: statistical and numerical approach. Int. J. Solids Struct. 40 (13–14), 3647–3679 . Kosmann, N. , Karsten, J. , Schuett, M. , Schulte, K. , Fiedler, B. , 2015. Determining the

effect of voids in GFRP on the damage behaviour under compression loading using acoustic emission. Compos. Part B 70, 184–188 .

Maragoni, L. , Carraro, P. , Peron, M. , Quaresimin, M. , 2017. Fatigue behaviour of glass/epoxy laminates in the presence of voids. Int. J. Fatigue 95, 18–28 . Maragoni, L. , Carraro, P. , Quaresimin, M. , 2016. Effect of voids on the crack formation

in a [45/45/0]s laminate under cyclic axial tension. Compos. Part A 91, 493–500 . CompTest 2015

Nguyen, V.P. , Lloberas-Valls, O. , Stroeven, M. , Sluys, L.J. , 2010. On the existence of representative volumes for softening quasi-brittle materials - a failure zone av- eraging scheme. Comput. Methods Appl. Mech. Eng. 199 (45–48), 3028–3038 . Park, K. , Paulino, G.H. , 2011. Cohesive zone models: a critical review of traction-sep-

aration relationships across fracture surfaces. Appl. Mech. Rev. 64 (6) . Saeb, S. , Steinmann, P. , Javili, A. , 2016. Aspects of computational homogenization at

finite deformations: a Unifying review from Reuss’ to Voigt’s bound. Appl. Mech. Rev. 68 (5) .

Scott, A. , Sinclair, I. , Spearing, S. , Mavrogordato, M. , Hepples, W. , 2014. Influence of voids on damage mechanisms in carbon/epoxy composites determined via high resolution computed tomography. Compos. Sci. Technol. 90, 147–153 . Suiker, A.S.J. , Turteltaub, S. , 2007a. Crystalline damage growth during martensitic

phase transformations. Philos. Mag. 87 (32), 5033–5063 .

Suiker, A.S.J. , Turteltaub, S. , 2007b. Numerical modelling of transformation-induced damage and plasticity in metals. Modell. Simul. Mater. Sci. Eng. 15 (1, SI), S147–S166 . IUTAM Symposium on Plasticity at the Micron Scale, Tech Univ Den- mark, Lyngby, DENMARK, MAY 21–25, 2006

Summerscales, J. , 1998. Microsctructural Characterisation of Fibre-Reinforced Com- posites. Woodhead Publishing Limited, Cambridge, England .

Svenning, E. , Fagerström, M. , Larsson, F. , 2016a. Computational homogenization of microfractured continua using weakly periodic boundary conditions. Comput. Methods Appl. Mech. Eng. 299, 1–21 .

Svenning, E. , Fagerstrom, M. , Larsson, F. , 2016b. On computational homogenization of microscale crack propagation. Int. J. Numer. Methods Eng. 108 (1), 76–90 . Turteltaub, S. , van Hoorn, N. , Westbroek, W. , Hirsch, C. , 2018. Multiscale analysis of

mixed-mode fracture and effective traction-separation relations for composite materials. J. Mech. Phys. Solids 117, 88–109 .

Unger, J.F. , 2013. An FE2-X1 approach for multiscale localization phenomena. J. Mech. Phys. Solids 61 (4), 928–948 .

Vajari, D.A. , 2015. A micromechanical study of porous composites under longitudi- nal shear and transverse normal loading. Compos. Struct. 125, 266–276 . Vajari, D.A. , Gonzalez, C. , Llorca, J. , Legarth, B.N. , 2014. A numerical study of the

influence of microvoids in the transverse mechanical response of unidirectional composites. Compos. Sci. Technol. 97, 46–54 .

Vajari, D.A. , Legarth, B.N. , Niordson, C.F. , 2013. Micromechanical modeling of uni- directional composites with uneven interfacial strengths. Eur. J. Mech. A 42, 241–250 .

Xia, Z. , Curtin, W. , Peters, P. , 2001. Multiscale modeling of failure in metal matrix composites. Acta Mater. 49 (2), 273–287 .