Fiber neutrality in fiber-reinforced composites

Fiber neutrality in fiber-reinforced composites

Evidence from a computational study

Goudarzi, M.; Simone, A.

DOI

10.1016/j.ijsolstr.2018.07.023

Publication date

2019

Document Version

Final published version

Published in

International Journal of Solids and Structures

Citation (APA)

Goudarzi, M., & Simone, A. (2019). Fiber neutrality in fiber-reinforced composites: Evidence from a

computational study. International Journal of Solids and Structures, 156–157, 14-28.

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

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

International

Journal

of

Solids

and

Structures

journalhomepage:www.elsevier.com/locate/ijsolstr

Fiber

neutrality

in

fiber-reinforced

composites:

Evidence

from

a

computational

study

M.

Goudarzi

a,∗

_,

_A.

_Simone

a,b

a Faculty of Civil Engineering and Geosciences, Delft University of Technology, Delft, the Netherlands b Department of Industrial Engineering, University of Padova, Padua, Italy

a

r

t

i

c

l

e

i

n

f

o

Article history:

Received 1 February 2018 Revised 29 May 2018 Available online 1 August 2018

Keywords:

Fiber-reinforced composites Embedded reinforcement Fiber neutrality

a

b

s

t

r

a

c

t

Wereportnumericalevidenceforneutralityofthinfiberstoaprescribeduniformstressfieldina

fiber-reinforcedcomposite.Elasticfiniteelementanalysesoffiber-reinforcedcompositesarecarriedoutwith

aconventionalfully-resolved modeland anoveldimensionally-reducedfiber model.Thetwomodeling

approachesarecomparedintheanalysisofmechanicalpropertiesandmatrix-fiberslipprofiles.An

anal-ysisoftheeffectivenessofvariousfiberorientationswithrespecttotheloadingdirectionshowsthatthe

notionofinclusionneutrality,originallyformulatedforrigidlineinclusionsbyWangetal.[Journalof

Ap-pliedMechanics,52(4),814–822,1985],holdsalsoforlinearelasticthinfiberswithimperfectinterface.

© 2018 The Authors. Published by Elsevier Ltd.

ThisisanopenaccessarticleundertheCCBY-NC-NDlicense.

(http://creativecommons.org/licenses/by-nc-nd/4.0/)

1. Introduction

Apart from fiber shape, surface treatment, and volume frac-tion, fiber spatial orientation is an important characteristic con-trollingload-transmissionmechanismsinfiber-reinforced compos-ites(KangandGao,2002).Althoughtheeffectoffiberorientation can be accurately assessed by means of computational homoge-nizationtechniques(Bergeretal.,2005;Xiaetal.,2003;Mortazavi etal.,2013;Shengetal.,2004;Lusti andGusev,2004),the gener-ationofaconformalfiniteelementmeshforcompositeswiththin fibersisa tediousandtime-consumingtask.Embedded reinforce-menttechniques,inwhichfiberdiscretizationisindependentfrom thediscretizationofthecompositedomain,canbeeffectivelyused forthisclass ofproblems.Here we assessthe validity ofa novel embeddedformulationandemployitto showinherent character-isticsofcompositesreinforcedwiththinfibers.

Generallyspeaking,twoclassesofmethodsareavailableforthe micro-mechanicalstudyoffiber-reinforcedcomposites:mean-field anddirectnumericalmethods.Althoughmean-fieldmethodssuch asEshelby-basedtwo-stephomogenizationschemes(Pierardetal., 2004; Tianetal., 2016) are fastand cost-effective,finiteelement (FE)averaging methods havegainedpopularity fortheir accurate

∗ _{Corresponding author.}

E-mail addresses: m.goudarzi@tudelft.nl (M. Goudarzi), angelo.simone@unipd.it ,

a.simone@tudelft.nl (A. Simone).

geometricalrepresentationofthecompositemicro-structure.A se-rious drawback of classical FE-basedhomogenization is the con-formal meshgenerationprocess forcompositeswithmanyfibers. When the number of fibers is relatively low, the composite can still bediscretizedusingclassical conformalapproachesasshown by Lusti and Gusev (2004) and Tian et al. (2016): in Lusti and Gusev (2004, Fig. 1) the discretization for a composite with350 nanotubes ofaspect ratio200 ata volume fractionof 0.5% con-sists of 3.5× 106 _{nodes tessellated into 21× 10}6 _tetrahedral

ele-ments;inTianetal.(2016),theauthorsshow thatthey can gen-erateadiscretizationforacompositefibervolumefractionsupto 20%considering600fibers(usingthereforefiberswitharelatively largediameter).Obviously,thestudyofcompositeswiththousands ortens of thousandsof fibers withvolume fractions around 20% (typical ofnanocomposites Andrews et al.(2002)) wouldbe un-feasible with conformal approaches. Advanced discretization ap-proachessuch asembeddedreinforcement methods(Balakrishnan andMurray, 1986; Elwi andHrudey, 1989; Hartl, 2002; Barzegar andMaddipudi,1997;Nini´cetal.,2014;Radtkeetal.,2010) facili-tatethediscretizationofhighaspectratiofibersbyallowingtheir mesh-independentrepresentation.Intheliterature,thesemethods have beenapplied to curved inclusions(Elwi andHrudey, 1989), slidingfibers(BarzegarandMaddipudi, 1997;Hartl,2002; Radtke etal.,2010), andtodescribetheinteraction betweenpile founda-tions andthesurrounding soil(Nini´c etal., 2014).In the

simula-https://doi.org/10.1016/j.ijsolstr.2018.07.023

tion in Section 5we have used up to around 22000 fibers with relativelycoarsediscretizations.

To lift the conformal meshing requirements of standard FE methods, two key assumptions are made in embedded formula-tions:(i)highaspectratiofibersaredescribedasone-dimensional barorbeamelements;and(ii)thefiberkinematics isintroduced by means ofa displacementgapbetween matrixandfiberwhile preserving the continuity of the underlying matrix displacement field across the fiber. Relaxing the second assumption in an em-bedded formulation requires the use of special enriched approx-imations (Pike and Oskay, 2016) that would significantly reduce the cost-effectiveness of the method and increase its complex-ity. For the analysis of composites with many fibers we prefer to adopt the most convenient approach from the computational point ofview.In Section2 we thereforeproposea novel embed-ded reinforced technique, which is compared to a fully-resolved three-dimensional fiber-reinforced model in Section 3. The nu-merical formulationproposed inthispaperissuperior toexisting embedded reinforcement models in that it uses totally indepen-dent finite element meshesfor fiberandmatrix. Thiswill elimi-natetheneedforcalculatingintersectionsbetweenfiberand ma-trix elements, which may require complex algorithms especially forcurvedfibers (Durandetal., 2015).Evenif thecomputational setup isverysimple,theone-to-onecomparisonbetweena fully-resolvedandadimensionally-reducedfibermodelallowsusto as-sesstheaccuracyofthedimensionally-reducedmodelintermsof global(homogenizedeffectivestiffness)andlocal(matrix-fiberslip profiles) quantities.Tothe author’sknowledge, thistype of com-parisonhasnotbeenreportedinpreviousstudiesandvalidatesthe use of dimensionally-reduced models for fibers with sufficiently smalldiameters.

Numerical homogenization studies of fiber orientation ef-fects(Tianetal.,2016;2014;KangandGao, 2002)show thatthe composite stiffness is very sensitive to fiber misalignments. In a numerical study, Tian et al. (2016, Fig. 7b) showed that the ef-fective Young’s modulus along the loading direction in a carbon fibercompositedecreasesratherquicklywhenthefibersare mis-alignedwithrespecttotheloadingdirectionandreachesits min-imum value, witha decreaseof ≈ 25%fromits maximum,when thefibersareinclinedat60° totheloadingdirection.Althoughthe authors didnot establisha correlation betweenelasticproperties andgeometricalpropertiessuchasfiberorientationanddiameter, their resultscan be related, mutatis mutandis, to those obtained for zero-thickness rigid inclusions (or rigid line inclusions, RLIs) problems(Wangetal., 1985).In RLIproblems,thesolution fields showastrongdependenceontheinclusionorientationandmatrix Poisson’s contraction effects,with the limit caseof the inclusion beingneutral(this,e.g.,meansthatthestressfieldisnotperturbed forcertaininclusionorientations).Intheliterature,inclusion neu-tralityhasbeendemonstrated onlyforperfectlybondedRLIs and for in-plane statesin experiments (Noselliet al., 2010), theoreti-cal studies(Wangetal.,1985;DalCorsoetal.,2016),and simula-tions (BarbieriandPugno, 2015). Inthis contributionwe demon-strate for the first time the validity of the concept of inclusion neutrality forhighaspect ratiofibers ina three-dimensional vol-ume. In Section 4,we show variousforms ofinclusion neutrality inafiber-reinforced compositewithlinearelasticthinfiberswith perfect and imperfect matrix-fiber interface using fully-resolved and dimensionally-reduced fiber models. Results are reported in terms of Young’s moduli, shear moduli and Poisson’s ratios. To demonstratetheapplicabilityofthedimensionally-reducedmodel tocompositeswithmanyfibers,Section5reportsadetailed micro-mechanical study. Results, including some observations on fiber neutralitywithrespecttotheeffectivePoisson’sratioofthe com-posite, arecompared withrule-of-mixturespredictions,whilethe limitationsofthelatterarepointedout.

2. Method

2.1. Weakformofthegoverningequations

Withreferencetotheprinciple ofvirtualwork forasmall dis-placementelastostatics problem, the weakform ofthe governing equationsforafiber-reinforcedcompositewithoutbodyforcesis  m

σ

m:

∇

s

δ

umd



m+  f

σ

f:

∇

s

δ

ufd



f+  int tc·

δ

wd



int −  t ¯t·

δ

umd



t=0, (1)

where

∇

s_{isthesymmetric-gradientoperator,andtheintegralover}

thetotalvolume



=



m∪



f issubdividedintomatrixandfiber

contributions. The virtual displacement vectors

δ

um and

δ

uf and

thecorresponding stresstensors

σ

m and

σ

f aredefinedover

ma-trix(m)and fiber(f)regions. The contributionalong the matrix-fiber interface



int represents the virtual mechanical work done

bythe interfacetractions tc forthelocalvirtual opening(sliding)

vector

δ

wacross(between) thetwosidesoftheinterface.The in-terfaceintegral isevaluated overthe interface surface



int of the

embeddedfibersand,forconvenience,isexpressedinacoordinate systemlocaltoeachfiber.Thelasttermin(1)istheworkdoneby theexternaltractions ¯tovertheexternalsurface



t.

2.2.Discretizedweakform

In the context of the finite element method, the domain



isdiscretized into matrixandfiber finiteelements, each withits ownsetofdegreesoffreedom.Tointroducefiberandmatrix-fiber interface contributions (second andthird terms in(1)), reference ismadeto two fiberdiscretizations:aconformal approachanda mesh-independentapproach.Whiletheapproximationofthefiber displacementvectoruf andtheopening/slidingvectorw,later

re-ferred to as the interface gap vector, is specific to each scheme, thefinal discretizedformofthegoverningequationsisvery sim-ilar for both techniques. Next, the discretization of matrix, fiber, andmatrix-fiberinterfaceisdiscussedindetail.

2.2.1. Matrix

The matrixdisplacement field um is approximated atthe

ele-mentlevelas um

(

x

)

≈ uhm

(

x

)

= n  i=1 N i

(

x

)

ui =Nud, (2)

whereNi(x) andui arethestandard Lagrangeshape functionand

displacementvectorsdefinedatnodei,respectively,thenumberof nodesinanelementisn,thematrixNu containselementalshape

functions,anddisthenodaldisplacementvector. Stiffnessmatrix Kme =

 e

m

BTuDmBud



m (3)

andexternalforcevector fexte =

 e

t

NuT¯td

,

(4)

withDmthematrixofelasticconstantsandButhematrixofshape

function derivatives, are obtained following standard procedures appliedtothefirstandlasttermsin(1).

2.2.2. Solidfibermodel

Areferencecomputationalmodelisdevelopedbymeshing the exactgeometryofthefiberin



fwithsolid-typeelementsandby

usinga conformal discretization atthe matrix-fiber interface. Al-thoughhighlyaccurate,thisapproachrequirestheconstructionof

aconformal discretization,atediousoperation formanypractical problems.Here,bothmatrixandfiberregionsarediscretizedusing conventionalnon-structuredlineartetrahedral elementsasshown inFig.1aforaninclinedcylindricalfiber.

Usingthecounterpartof(2)forthefiberdisplacementfield,the stiffnesscontributionrelatedtothesecondtermin(1)is

Ke f =  e f BT uDfBud



f, (5)

whereDf isthematrixofelasticconstantsofthefibermaterial.

Todiscretizetheinterfacegapvectorw,zero-thickness confor-malinterfaceelementsareincludedbetweenmatrixandfiber.This requiresfindingallmatrixelementssharingatleastonenodewith fiberelements(shaded in red inFig. 1a); shared nodesare then duplicated and zero-thickness interface elements are embedded betweenmatrix and fiberas depicted inFig. 1b. The discretized interfacegapvector,

wh_=RN

intdint, (6)

isexpressedinacoordinatesystemlocaltotheinterface.Without lossof generality and by making referenceto the node number-ing inFig. 1b,this is achievedby defining the rotation matrix R

fromtheglobaltothelocalcoordinatesystem,theinterfaceshape functionsmatrix Nint=



Nf int

|

N m int



=



N 1I N 2I N 3I

|

−N4I −N5I −N6I



(7) expressed in terms of isoparametric coordinates for the zero-thicknessinterface element,andtheinterface nodaldisplacement vector dint=



uf1 uf2 uf3

|

um4 um5 um6



T , (8)

where I is the identity matrix of size three, and u_fi(i=1,2,3) and

umi(i=4,5,6)indicatethedisplacementvectorsforfiberandmatrix, respectively,intheglobalcoordinatesystem.

Followingstandardprocedures,theinterfacestiffness contribu-tion Ke int=  int NT

intRTDbRNintd



int (9)

isobtainedonsubstituting(6)intothethirdtermin(1)andwith the interface constitutive matrix Db= ∂_∂twc =diag

(

Kbt,Kbn1,Kbn2

)

fora linear elastic interface. The constantsKbn1 andKbn2

repre-sentthestiffnessesoftheinterface inthedirectionnormaltothe interfacesurfaceandperpendicular toitandtothefiberaxis, re-spectively,andKbtrepresents thestiffnessof theinterface inthe

directiontangentialto thefiberaxis.In thiswork,we onlyallow fiberslipinthedirectiontangentialtothefiberaxisby constrain-ingthenormaldisplacementgapsbymeansofpenalty augmenta-tionofthe interface stiffnessesKbn1 andKbn2.Additionally,

inter-faceelementsarenotintroducedatthefiberendpointswherewe assumenoadhesionbetweenfiberandmatrix.

Expansionoftheinterfacestiffnesscontribution(9),considering thattheinterface gapvector (6)isa functionofboth matrixand fiberdisplacements,resultsin

Ke int=



Kmm int K mf int Kfm int K ff int



=



 intN m int T RT_D bRNintmd



int  intN m int T RT_D bRNintf d



int  intN f int T RT_D bRNintmd



int  intN f int T RT_D bRNintf d



int



. (10)

2.2.3. Dimensionally-reducedfibermodel

When modeling highaspectratio fibers,theconforming finite elementmodeldescribedinthe previoussection can bereplaced byamesh-independentdimensionally-reducedmodelwitha dras-ticreductionofthecomputationalcost.

In the proposed model, inspired by embedded reinforcement techniques (Balakrishnan and Murray, 1986; Elwi and Hrudey, 1989;BarzegarandMaddipudi,1997),fibersareidealizedas one-dimensional objects that respondto axial deformation only. This is a reasonable assumption when dealing with high aspect ra-tiofibers.Afiberdiscretizationindependent fromthematrix dis-cretizationandconsistingoflinearone-dimensionalLagrange ele-mentsisusedtorepresentembeddedfibers.Fig.2showsthecase ofasingleembeddedfiberalthoughtherearenolimitationsonthe numberoffibersthatcancrossa matrixelement.Worthnoticing istheindependenceofthetwodiscretizations:fibernodesdonot coincide withintersectionpointsbetweenfiberaxis andelement faces.Whileanystructuredorunstructuredfiberdiscretizationcan beused,inthenumericalsimulationsinthispaperfibernodesare uniformlydistributedforconvenience.

Tosimplifythegenerationofthe stiffnessmatrix forelements thatareintersectedbyone-dimensionalfibers,stiffnessmatrix(3), with the integration performedover the total volume



e _{of the}

element,isusedforthematrixcontribution.Thishoweverrequires theuseofaneffectiveYoung’smodulusforthefiberstocancelout thealreadyconsideredmatrixcontributioninthefiberregion.The contributionofthekthfiberelement(withendpoints aandbfor theredfiberelementinFig.2) isidenticalto thestiffnessmatrix ofaone-dimensionalstandardtrusselementinthree-dimensional spaceandiswrittenas

Kk f =

(

E k f − Em

)

A k L k



Tk _−Tk −Tk _Tk



, (11)

whereEmandE_fkarethematrixandfiberYoung’smoduli,

respec-tively,Ak_{is thecrosssectionalarea(A}k₌

_π

_d2

f/4withdfthe fiber

diameter),andLk_{isthefiberelementlength.Thefiberorientation}

matrix Tk₌

⎡

⎣

e k xe kx e kxe ky e kxe kz e k xe ky e kye ky e kye kz e k xe kz e kye kz e kze kz

⎤

⎦

(12) with e k x= x k b− x k a L k , e k y= y k b− y k a L k , e k z= z k b− z k a L k , (13)

where (xa,ya, za) and(xb, yb,zb) indicate the coordinates ofthe

fiberelementends(e.g.,nodesaandbfortheredfiberelementin Fig.2).

Consequently,thematrix-fiberinterfaceisreducedtoan equiv-alentlineobject,withinterfacegapvector(6)andstiffnessmatrix contribution(9)evaluatedalongit.The interfacegapvector(6)is definedwith Nint=



Nf int

|

Nintm



=



N f aI N bf I

|

−N1m I · · · −NnmI



(14) and dint=



uf a ufb

|

um1 · · · umn



T , (15) where Nf

a andNfb are the one-dimensional Lagrange shape

func-tionsatnodesaandbdefinedinanisoparametriccoordinate sys-temlocaltothefiber,andtheshapefunctionNm

i atnodei ofthe

parentmatrixelement isalsodefinedinanisoparametric coordi-natesystem.

Inasmalldeformationsetting,theinterfacegapvectoris eval-uatedintheundeformedconfigurationatintegrationpointswithin fiber elements (e.g., at a and b for the red fiber element in Fig. 2 since we consider a two-point Gauss–Lobatto quadrature scheme). Thisimplies that once thelocal coordinates ofan inte-grationpointalongafiberelementaredefined,theyareexpressed into theglobalcoordinate systemandused to identifyits parent matrixelementbymeansofanoctreestructure(Duflot,2006, Sec-tion9.1)thatisconstructedpriortotheanalysis.

The fiber shape function in (14) are directly evaluated at the integration point in an isoparametric coordinate system local to thefiber.Theshapefunctionofthematrixelementarealso evalu-atedatthesamepointnowexpressedintheisoparametric coordi-nate systemrelated tothe matrix elementand obtainedthrough an inverse mapping procedure starting from its position in the globalcoordinatesystem.Matrixandfiberelementsareconnected throughso-calledbondelementsthataredefinedwithreferenceto thetwofiberelementnodesandthenodesofthematrixelement inwhichthefiberelementnodesarelocated(thediscretizationis thereforeconformingontheside ofthe fiberandnotconforming on thesideofthe matrix,anda fiberelementcan havenodesin twodistinctmatrixelements).

The interface contribution(9)relatedto thekthfiberelement becomes Kk int=C k f  Lk NT intRTDbRNintds, (16) withCk

f =

π

dfthecircumferenceofthefiber,Lkthelengthof

fiber-matrix bond element, and sthe local coordinate along the bond element.Thus,forasystemwithn_f fiberseach subdividedintons

lineelements,asshowninFig.2,thetotalnumberoffiberand in-terfaceelementcontributionsisequalto2× nf× ns.Similarto(10),

theinterface contributionissharedbetweenmatrixandfiber de-grees of freedom in the assembly of the corresponding stiffness matrix.

2.2.4. Globalsystemofequations

Thenf× ns fiberelement contributionsandnf× ns matrix-fiber

interface element contributions to theglobal stiffnessmatrix are assembledseparatelyfromthoseofthematrixelements. Irrespec-tive ofthefiberdiscretizationscheme,the generalformofglobal systemofequationsis

⎡

⎢

⎢

⎢

⎢

⎣

Kmm Kmf1 Kmf2 . . . Kmfn_f Kf1m Kf1f1 0 . . . 0 Kf2m 0 Kf2f2 . . . 0 . . . ... ... ... ... Kfn_fm 0 0 . . . Kfn_ffn_f

⎤

⎥

⎥

⎥

⎥

⎦

⎡

⎢

⎢

⎢

⎢

⎣

m f1 f2 . . . fnf

⎤

⎥

⎥

⎥

⎥

⎦

=

⎡

⎢

⎢

⎢

⎢

⎣

Fm 0 0 . . . 0

⎤

⎥

⎥

⎥

⎥

⎦

, (17)

wheremandfi

(

i=1,...,nf

)

arevectorscontainingmatrixand

in-dividual fiberdisplacements, respectively. The submatricesin the globalstiffnessmatrixin(17)areassembledbyapplyingstandard procedurestoelementalstiffnesscontributionsofmatrix,fiberand matrix-fiber interface. In this framework, matrix and fiber stiff-nessesleadtoblockdiagonaltermswhile,asshownin(10), inter-facecontributionsensuringthecouplingbetweenmatrixandfibers leadtobothblockdiagonalandblockoff-diagonalterms.

Throughoutthisstudy,linearelasticbehaviorisassumedforthe constitutivemodelsofmatrix,fiber,andmatrix-fiberinterface.

3. Validityofthedimensionally-reducedapproach

Theinfluenceoftheone-dimensionalreductionofthefiber ge-ometry onthemechanicalresponseofa fiber-reinforced compos-ite is studied. The geometrical reduction is performed over the fiberdiameter df,with fixedfiber lengthlf. The validationofthe

dimensionally-reduced fibermodelis performedagainstthe solid fibermodeldescribedinSection2.2.2assumingtheratiobetween fiberdiameteranddomainsize Lasthecharacterizingparameter. To reduce thecomplexityof thediscretization procedures dueto thegenerationofaconformingmesh,onlyonefiberisconsidered intheconfigurationshowninFig.1aforthereferencemodel. Un-lessstatedotherwise,inall thesimulations inthispaperwe have assumedthat matrix Young’s modulus Em andPoisson’sratio

ν

m

are equal to 100 MPaand 0.4, respectively,and that inthe solid

Fig. 1. (a) Reference conformal finite element model with a cylindrical fiber that is discretized using tetrahedral finite elements. (b) Zero-thickness conformal inter- face elements are placed between fiber (blue region) and surrounding matrix (red region) to allow fiber slip. It is assumed that no adhesion exists between the end- points of the fiber and the matrix. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

fibermodel, matrixandfiber havethe same Poisson’sratio. Fur-ther,theeffectivemechanicalpropertiesofthecompositeare esti-matedthroughtheproceduredescribedinAppendixA.

3.1. Perfectbond:Effectivemechanicalproperties

We now determine the range of fiber diameters over which thepredictionsofthedimensionally-reducedfibermodelarevalid. Tothis end, we compute the values of the effective longitudinal Young’s modulus Ec

x of a composite with one fiber aligned with

thex-axis(

θ

z=0◦inFig.3)andcomparethemtothe

correspond-ingvaluesobtainedwiththesolidfibermodel.Itshouldbe men-tionedthatthisorientationischosenamongallorientationsinthe xyplanebecauseityieldsthelargestincreaseinEc

x.Computations

aredonewithtwodifferentfiberlengths(l_f=0_.4Landl_f=0_.8L).A largevalueoftheinterfacetangentialstiffnessKbtisusedtomimic

aperfectlybondedinterface(hereafterreferredtoasperfect inter-face).Fig.4showstheeffectivecompositeYoung’smodulusEc

x for

variousvaluesofthefiberYoung’smodulusE_f,bothnormalizedby thematrixYoung’smodulusEm.Byreducingtheratiodf/L,the

re-sponsesofthenumericalmodels convergetothe samevalue.For d_f/L_<0.01averygoodagreementbetweenthemholds,anda one-dimensionalrepresentationofthefibercanbejustified.

3.2.Imperfectbond:Matrix-fiberslip

Next,we studyhowtheslipprofilesbetweenmatrixandfiber compare forthe reduced modeland the solid fiber model when a linear elastic traction separation law is employed. While the dimensionally-reduced fiber model obviously generates one slip profile, an infinite number of slip profiles can be sampled with thesolid-fibermodel.Toobtainavisually meaningful representa-tionoftheresults,wehaveconsideredaveryweakinterface(K_bt= 500N/mm3),hereafterreferredtoasimperfectinterface.The sim-ulationboxisdeformedinthehorizontaldirectionthroughthe im-positionofaperiodicstrain(

ε

xx=0.1)asdiscussedinAppendixA.

Thiscomputational setup eliminates the influenceof the domain size Lfrom the results.Fig. 5a (5b) showsthe resulting normal-izedfiber slipprofiles fora fiber alignedwiththe loading direc-tion (

θ

z=0◦) and fiber diameterdf=0.03L (0.005L). With

refer-encetothe solid fibermodel,inthisspecial casethepoints sur-rounding the fiberat the same horizontalcoordinate experience the same displacement with respect to the corresponding fiber points.Auniqueslipprofileisthereforegeneratedwhichisfound to be ina remarkablygood agreement withthe slipprofile gen-erated by the dimensionally-reduced fiber model. In contrast, as shownin Figs. 5cand d,an inclinedfiber (

θ

z≈ 28.6◦

(

=0.5rad

)

)

Fig. 2. A dimensionally-reduced fiber intersects a solid matrix element: (a) three-dimensional view, (b) top view. The fiber can be placed within the solid element region and can be discretized independently from it. Discretized fiber elements are shown in blue, separated by nodes. The figure shows also a fiber element, in red, between nodes

a and b . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 3. An L × L × L ( L = 1 mm) periodic simulation box with a fiber of length l f in

the middle aligned along the x axis. In the fiber orientation study ( Section 4.1 ), the fiber can rotate around the z axis by an angle θz in the xy plane.

doesnot slide uniformly with respect to the surrounding mate-rial. The gray shaded region represents all the fiber slip profiles sampledover thelateralsurfaceindirectionsparallelto thefiber axis.Thesedisplacement slipprofiles cannot be predictedby the reduced model where, instead, a single slip profile is predicted (dashedlines),whichishoweveringoodagreementwiththe av-erageslipprofileofthesolidfibermodel.

As shown in Figs. 5c and d, the width of the shaded region decreases by reducing the fiber diameterindicating that the dif-ference betweentheupperandlower interface displacement gap vectors becomes smaller. This means that the results obtained withthesolidfibermodelconvergetothoseofthe dimensionally-reducedfibermodelasthefiberdiameterdecreases.Togeneralize the relation betweenfiber diameter and interfacialdisplacement gapsforinclinedfibers,thenormalizedarea Ae oftheshaded

re-gion,referredtoasnormalizedgaparea,isplottedagainstthe nor-malizedfiberdiameterin Fig.6. Resultsare shownin Fig.6a for differentvaluesofthefiberYoung’smodulusEfandfiberlengthlf.

Irrespectiveoffiberlengthandmaterialproperties,thenormalized gapareadecreaseswithdecreasingfiberdiameter,anditsvalue is obviouslyafunction offiberorientation,decreasingwith decreas-ingfiberorientationangleasillustratedinFig.6b.

4. Fiberneutrality

According to previous studies on planar reinforcements, perfectly bonded rigid inclusions become mechanically neu-tralunderspecific circumstances—mechanically neutralinclusions do not influence the mechanical response of the composite. Wangetal.(1985)havederivedtheanalyticalsolutionfora zero-thickness RLI subjected to an inclined loading at infinity; for a

given in-plane problem (plane stress/plane strain), the angle be-tween the loading direction and the RLI at which the inclusion does not influence the stress field depends only on the Pois-son’s ratioof thematrix material.A similarproperty hasnot yet beenreportedinfiber-reinforcedcompositesdespitevarious stud-iesontheroleoffiberorientationintheir homogenized mechan-icalproperties(Tianetal.,2016;2014;KangandGao,2002). Mo-tivated bytheseobservations, we performa studyonorientation effectsandstress neutralitypropertiesforlinearelasticthin(high aspectratio)fiberswithimperfectmatrix-fiberinterface.

4.1. Neutralityinthedimensionally-reducedmodel

Neutral orientations are determined for effective Young’s and shearmoduliandPoisson’sratios.Forthesakeofillustration, spec-imens withfiber volume fraction

υ

_f=1_.0% (1300 aligned fibers) areconsidered.Lengthanddiameterofthefibers,discretizedwith the dimensionally-reduced model,are set to 0.2L and0.007L, re-spectively.Theresultshavebeenobtainedwiththenumerical ho-mogenization procedure described in Appendix A. Numerical ef-fective properties of the homogenizedcomposite are reported in Fig.7aforfiberorientations

θ

zrangingfrom0° to90°.Inthe

sim-ulations, the ratioE_f/Em is set to 100 andan imperfect interface

betweenfiberandmatrixisassumed.Althoughwepresentresults forthisspecific setofmaterial properties,their validityforother parameter sets, also including the case of a rigid interface, was confirmedwithsimulationswhoseresultsarenotreportedhere.A comparisonwithanalyticalmicromechanicalestimatesisprovided inSection5.

Asexpected, thehighesteffectiveYoung’smodulusisobtained when fibers are aligned with the loading direction. Contrary to thegeneralunderstandingthatfibersalwaysincrease thestiffness of a fiber-reinforced composite, corroborated also by results ob-tainedwithtwo-stepmean-fieldhomogenizationprocedures (see, e.g.,Pierardetal.(2004,Fig.3)andTianetal.(2016,[Figs.7–10]), the minimum composite stiffnesscorresponds to a state of fiber neutrality.ItisadditionallyimpliedfromFig.7athatfibers perpen-diculartotheloadingdirection(

θ

z=90◦)can,ingeneral,increase

theeffectiveYoung’smodulusofthecomposite.

Fiberorientationscorrespondingtoaneutralsituation(i.e., neu-tral orientations

θ

n at which fibers have no influence on the

Young’s modulus ofthe composite) are function of the Poisson’s ratio

ν

m ofthe matrix material only andare reported inFig. 7b

(blue-filled circles).Inthesamefigure wereport valuesforplane states corresponding to the cases of unperturbed stress fields aroundazero-thicknessperfectly-bondedrigidplanarinclusion ac-cordingtotheanalyticalsolutionbyWangetal.(1985).Itis

inter-estingtonoticethatdimensionally-reducedfibersaremechanically neutralatthe neutralorientationsvalidforRLIsin aplane stress state.

As shown in Fig. 8, neutrality is observed in terms of the effective shear modulus Gc

xy of the composite as well. In this

case, neutral orientations are 0° and 90°. A similar stress neu-tralityhas beenreportedin RLIsunder uniformshearloading by Nosellietal.(2010)usingphotoelasticityinatwo-dimensional ex-perimentalstudy.Finally,variouscasesofneutralityintermsofthe effectiveanisotropicPoisson’sratioscanbeidentifiedinFig.9.

These results confirm that linear elastic thin fibers share the neutrality feature ofzero-thickness RLIs (i.e., theneutrality angle dependsonlyonthematrixPoisson’sratio).Asalreadymentioned, neutrality holds irrespective of fiber-matrix interface parameter. This is atvariance with previously known cases of neutrality of arbitraryshapedinhomogeneities(Ru,1998;BenvenisteandMiloh, 1999)inwhichneutralityisachievedbydesigningnon-idealor im-perfectinterfacesbetweeninclusionandmatrix.

4.1.1. Fibercompressionundertensileloading

Atneutralorientations

θ

nshowninFig.7b,afiberdoesnot

ex-perience deformation. Thisimpliesa signshift inthe fiberstrain andslip. Fig. 10showsthe slipandstrain profiles inthefiber of the one-fiber composite depicted in Fig. 3 under tensile loading along the x-axis at various orientations (neutral angle

θ

n≈ 55°).

Results are shownforthecaseof anearly incompressiblematrix with

ν

m=0.49; thevaluesof allother propertiesare takenfrom

theexampleintheprevioussection.Asimilarsetofresultscanbe obtainedforothermatrixPoisson’sratios.

Ashiftinsignwhen

θ

z exceedstheneutrality angle

θ

n canbe

detectedinbothfiberslipsandstrainprofiles.Tensilestrainsoccur infiberswithorientationsbelow

θ

nwhilefororientations

exceed-ing

θ

n the sign reverses andcompressive strains along the fiber

axisaregenerated.Thisbehaviorcouldtriggercompositefailurein theformoffibermicro-buckling(BudianskyandFleck,1993),even withexternallyappliedtensileforces.

4.2. Neutralityandfiberdiameterinthesolidfibermodelwith perfectinterface

In the previous section we have shown, using the dimensionally-reduced fiber model, that linear elastic thin fibers with imperfect matrix-fiber interface are neutral under specific circumstances. We now investigate how the fiber diameter in-fluences thisproperty considering a perfectlybonded solid fiber. Thecaseofan imperfectmatrix-fiberinterface isdiscussedinthe next section. The results are obtained using the simulation box inFig.3 withone perfectlybondedfiber(lengthl_f=0_.8L) that is discretized using the solid and the dimensionally-reduced fiber models.

Toinvestigatetherelationbetweenfiberdiameterandthe oc-currence of fiber neutrality, the fiber is oriented along the neu-tralorientationthatwasdeterminedforthedimensionally-reduced fibermodel(

θ

z=

θ

n≈ 57◦fromFig.7b,foramatrixwithPoisson’s

ratio

ν

m=0.4).ThenormalizedeffectiveYoung’smodulusExcis

re-portedinFig.11atdifferentfiberdiametersandstiffnesses.Forthe sakeofcomparisonwiththelongitudinalYoung’smodulus,which givesthemaximumachievableeffectivestiffness,solidfibermodel resultsforthecasewith

θ

z=0◦ arealsoincludedintheplots.

AsshowninFig.11,thedimensionally-reducedmodelpredicts neutrality at all fiber diameters when

θ

z≈ 57° (black curves in

Fig. 11). Since the fiber is described as a line object, this is at-tributed to the fact that the diameteris not geometrically incor-poratedinan explicitmannerbutisconsideredasaparameterin thefiberandinterfacecontributionstothestiffnessmatrix. How-ever, whenthe exactgeometry ofthe fiberistakeninto account

usingthesolidfibermodel(bluecurvesinFig.11),neutralityholds onlyforrelativelysmallfiberdiameters(df/L<0.01).Fordf/L>0.01

theincreaseoftheeffectiveYoung’smodulusismodestandisnot significantly influenced by the fiber stiffness unlike the effective longitudinalYoung’smodulus(redcurves).

Analogous observations hold for the effective shear modu-lus Gc

xy. Fig. 12 shows a comparison of the two fiber models at

θ

z=0◦ (according to Fig. 8, neutrality in terms of Gcxy occurs at

θ

z=0◦). Results are now accompanied by the solid fiber model

predictionsfora fiberorientation (

θ

z=45◦) that yieldsthe

max-imumeffectiveshearstiffness.

These observations are somehow in agreement with previ-ousnumericalstudies offiberorientationinshortfiber-reinforced metal matrixcomposites (KangandGao, 2002;Tianetal., 2016). In these studies a modest increase of the effective Young’s and shear moduli was found for fibers inclined at

θ

z=60◦ and

θ

z=

0◦ in short carbon fiber composites (fiber with aspect ratio 15, df≈ 0.03L).Incontrast,intheresultsshowninFigs.11and12,

neu-tralityoffibersisexpectedatsmallerdiameters,andtheeffective stiffnessofthecompositewillnotchangecomparedtothematrix stiffnesscomputedattheorientationsshowninFig.7b.

4.3.Neutralityandimperfectinterface

Neutrality was already discussed in Section 4.1 under imper-fectinterfaceconditionswiththeanalysisrestrictedtothe single-interfaceassumptionofthedimensionally-reduced fibermodel.In thissection,wediscusshowneutralitycanchangeconsideringthe actualinterfaceinathree-dimensionalcompositewithasolidfiber andatwo-dimensionalplanarcompositewithathininclusion.

4.3.1. Solidfibersinathree-dimensionalcomposite

Animperfect interface leadstopartial load transfer,with dis-placement jumps occurring between fiber and matrix.According to thenon-uniqueness ofthe slip profile foran inclined fiberas discussed in Section 3.2, displacement differences between fiber edges(e.g.,thetopandbottomedgesshowninFig.5candd)can occur. Forrelatively smallfiberdiameters,however, displacement differencesbetweenedges becomesmallandtheslipprofilescan be approximated by a unique curve. Here we further assess this approximationandthesignificanceofthedifferentslipprofilesby emphasizingtheireffectonthedevelopmentofstressneutrality.

Similar toSection 4.2, acomposite withone fiber orientedat the neutral orientation (

θ

z=

θ

n≈ 57◦ from Fig. 7b,

ν

m=0.4) is

discretized using the solid fiber model. Values of Young’s mod-ulus Ec

x are plottedas a function of fiberdiameter in Fig.13 for

bothperfectandimperfectinterfaces.Forthesakeofcomparison, the Young’s moduli obtained with the solid fiber model with

θ

z=0◦(blacklines)arealsoincludedintheplotforbothinterface

conditions.Comparedtotheperfectinterfacecase(solidlines),an imperfectinterface(dashedlines)leadsingeneraltolowervalues oftheeffectivelongitudinalYoung’smodulusEc

x.Morespecifically,

whentheangle

θ

z≈ 57°,thereductionisonlypronounced at

rel-ativelylarge fiberdiameters;atsmallfiberdiameters(df/L<0.01)

the longitudinal Young’s modulus of the composite converges to thevalueoftheYoung’smodulusofthematrixand,withaperfect orimperfect interface, a neutral state ispredicted. When

θ

z=0◦

areducedstiffnessisdetectableevenatsmallfiberdiameters(see insetinFig.13).

To summarize, irrespective of the interface type, a fiber with asmalldiametercomparedtothecharacteristicdimension ofthe domainisnotstrainedwheninclinedat

θ

z=

θ

n≈ 57◦(neutral

ori-entationfor

ν

m=0.4)andhasthereforeanulleffectonthe

Fig. 4. Normalized effective com posite Young’s modulus as a function of fiber diameter for a fiber with aspect ratio between 4 and 266. Predictions of both models agree for small fiber diameters relative to the domain size.

Fig. 5. Normalized fiber slip s for a single embedded fiber with imperfect interface with E f / E m = 10 and l f = 0 . 8 L . (a,b) A unique slip profile can be identified for a fiber

aligned with the loading direction. (c,d) The slip profile is not unique for an inclined fiber. The dimensionally-reduced model produces a unique slip profile that agrees with the average slip of the solid fiber model (dashed line). Refer to Fig. 1 a for the nomenclature used for the solid fiber model.

4.3.2. Dimensionally-reducedinclusioninatwo-dimensional composite

The two-dimensional conformal model in Fig. 14,obtained as alimitcaseofthepreviouslydiscussedsolidfibermodel,can ade-quatelymodelcertaintypesoftwo-dimensionalcomposites(Sheng etal.,2004;Halletal.,2008).Themodelisemployedtoinvestigate theoccurrenceofneutralstates(reportedinWangetal.(1985)for perfectly-bondedtwo-dimensionalplanarinclusions)under imper-fect interface conditions. In this computational setting, two

in-terface surfacescan be identified when the inclusion is not per-fectly bonded to the matrix: one at each side of the inclusion and physically disconnected from each other. Depending on the out ofplane shapeofthe inclusion,differentassumptions can be maderegardingthecontinuityofthedisplacementfield acrossit. For a situation with platelet inclusions, as in clay nanocompos-ites (Sheng et al., 2004), displacement jumps across the platelet adequately represent the expected kinematics even for platelets of smallthickness. As soon asthe inclusion out-of-plane

dimen-Fig. 6. Effect of fiber diameter on the normalized gap area for a fiber with aspect ratio between 2 and 160. (a) Normalized gap area for fibers with different properties and

θz ≈ 28.6 ° ( = 0 . 5 rad). (b) Normalized gap area for different fiber rotationsθz .

Fig. 7. (a) Normalized effective Young’s modulus E c

x calculated at various fiber orientations θz and matrix Poisson’s ratios νm . (b) The neutral orientation θn is the orientation

corresponding to fibers experiencing no deformation under the applied load. Results obtained with E f /E m = 100 , K bt = 500 N / mm 3 , and υf = 1% (1300 aligned fibers with

df = 0 . 007 L and l f = 0 . 2 L ). Plane stress and plane strain results are from Wang et al. (1985) . (For interpretation of the references to colour in this figure legend, the reader is

referred to the web version of this article.)

Fig. 8. Normalized effective shear modulus G xy calculated at various fiber ori-

entations θz and matrix Poisson’s ratios νm . Results obtained with E f /E m = 100 ,

Kbt = 500 N / mm 3 , and υf = 1% (1300 aligned fibers with d f = 0 . 007 L and l f = 0 . 2 L ). sion decreases and the inclusion represents a fiber, as in fiber networks in ordinary paper sheets or buckypaper (Hall et al., 2008), the double-interface modelis not representative anymore, and a traditional single-interface model ofthe type describedin Section2.2.3istobepreferred.

In the double-interface model shown in Fig. 14, conceptually similartothemodelbyPikeetal.(2015),zero-thicknessinterface elementsareplacedateachsideoftheinclusiontoallowthe oc-currence of relative displacements betweeninclusion andmatrix on bothsides oftheinclusion. Inthis situation,upperandlower interfacesmoverelativetoeachother.Todescribethiskinematics, theinterfaceshapefunctionmatrix(7)isrewrittenas

N_intt/b=



N i 1 0 N i2 0 −N t/b 1 0 −N t/b 2 0 0 N i 1 0 N i2 0 −N t/b 1 0 −N t/b 2



, (18) where Ni

1 and N2i are one-dimensional Lagrange shape functions

evaluated atthe inclusion end points1 and2, respectively. Sim-ilarly, N₁t/b and N₂t/b are the matrix nodal shape functions evalu-atedat points1 and2,and superscripts t andb denotetop and bottominterfaces,respectively.Thesimplersingle-interfacemodel, whichdoes not account fora discontinuousmatrix displacement field,canbe recoveredby constrainingthe topandbottomnodes toexperiencethesamedisplacement(ut_=ub_).

Figs. 15a–c show the effective mechanical properties of the 1mm× 1mm periodic plate with a centered 0.5 mm long in-clinedinclusion asa functionoftheinclusionorientation

θ

z.

Me-chanical properties are extracted using the procedure described inAppendixA.Intheanalysesweusedtwovaluesofthe interfa-cialtangentialstiffnessK_bt(500N/mm2_and1000N/mm2_),andwe

assumedaunit thicknessfortheplateintheout-of-plane dimen-sion.TherelativelylowvaluesshowninFig.15areaconsequence ofthechoiceoftheparametersandthemodestmechanicaleffect ofasinglefiber.

Fig.15ashowsthatinthedouble-interfacemodeltheangle

θ

z

thatcorrespondstotheminimumeffectiveYoung’smodulusofthe composite (marked in the figure) is a function of the interfacial tangentialstiffnessKbt: thevalue ofthe angleincreases with

in-creasingvaluesofKbtandtendsto ≈ 50.7°.Forthesingle-interface

model,thevalue ofK_bthas noeffectontheangleatwhich Ec

x is

ataminimum: thisangleis equalto ≈ 50.7° andcorresponds to theneutralorientation previously discussed (inthis specificcase,

θ

n≈ 50.7° for

ν

m=0.4 under plane strain asshown in Fig. 7b).

This means that in the range of validity of the dimensionally-reducedmodel,Ec

x willassume itsminimumvalue,corresponding

toasituationoffiberneutrality,atthesameinclusionangleshown inFig.7b,irrespectiveofthevalueoftheinterfacialtangential stiff-nessKbt.

Fig. 9. Normalized effective Poisson’s ratio calculated at various fiber orientations θz and matrix Poisson’s ratios νm . Results obtained with E f /E m = 100 , K bt = 500 N / mm 3 ,

and υf = 3% (1300 aligned fibers with d f = 0 . 007 L and l f = 0 . 2 L ).

Fig. 10. A one-fiber composite with a nearly incompressible matrix material ( νm = 0 . 49 ) under the action of a horizontal external load. Fiber slip s (a) and fiber axial

strains εf (b) are strongly influenced by the fiber orientation θz . Results obtained with E f /E m = 100 , K bt = 500 N / mm 3 , d f = 0 . 007 L and l f = 0 . 8 L .

Fig. 11. Normalized effective Young’s modulus of the one-fiber composite (with l f = 0 . 8 L ) as a function of fiber diameter for solid and dimensionally-reduced models at

neutral orientation ( Fig. 7 b: θz ≈ 57 ° for νm = 0 . 4 ) of the dimensionally-reduced model for a fiber with aspect ratio between 8 and 266. The longitudinal Young’s modulus

of a composite with a horizontal fiber ( θz = 0 ◦) obtained with the solid fiber model (previously reported in Fig. 4 , second row) is also included. (For interpretation of the

Fig. 12. Normalized effective shear modulus of the one-fiber composite (with l f = 0 . 8 L ) as a function of fiber diameter for solid and dimensionally-reduced models at

neutral orientation ( θz = 0 ◦from Fig. 8 , irrespective of νm ) of the dimensionally-reduced model for a fiber with aspect ratio between 8 and 266. The maximum shear

modulus obtained with the solid fiber model for a fiber inclined at θz = 45 ◦is also included.

Fig. 13. Normalized Effective Young’s modulus of a composite with one fiber ori- ented at the neutral orientation obtained with the dimensionally-reduced fiber model ( θz ≈ 57 ° from Fig. 7 b, νm = 0 . 4 ) and discretized using the solid fiber model

under perfect and imperfect interface conditions for a fiber with aspect ratio be- tween 8 and 266. The longitudinal Young’s moduli obtained with the solid fiber model (previously reported in Fig. 4 , second row and second column) are also in- cluded. Both perfect and imperfect interface conditions are under a neutral state at small fiber diameters with θz = θn ≈ 57 ◦. Results obtained with E f /E m = 100 ,

Kbt = 500 N / mm 3 and l f = 0 . 8 L .

Fig. 14. Single conforming inclusion with elastic imperfect interface. Interface ele- ments are placed at each side of the inclusion. For illustrative purposes, both nor- mal and tangential displacement jumps are depicted.

In the single-interface model the angle corresponding to the minimumvalue ofEc

x indicatesa situationofinclusion neutrality.

Inthedouble-interfacemodelhoweverweobservethati)the com-positeYoung’smoduluscanbelowerthanthatofthematrixfora wide rangeofinclusion angles (stiffnessdegradation), andii)the amplitudeoftherangedecreaseswithincreasingK_btvalues. Simi-larconsiderationscanbedrawnfortheshearmodulusinFig.15b althoughneutralityisonlyobservedforthesingle-interfacemodel at

θ

equalto0° or90°.

Undershearloading,relativemovementsbetweenedges ofthe inclusion occur and, asshownin Fig. 15b,a similar reduction in thevalueoftheeffectiveshearmodulustoavaluelowerthanthe matrixshearmodulusisobtainedwiththedouble-interfacemodel. Additionally, maximum shear modulus is achieved for

θ

z=45◦,

wherepredictions ofsingle- anddouble-interface models are the same.

AsforthePoisson’sratioinFig.15c,thedouble-interfacemodel yields an increase of the composite effective Poisson’s ratios

ν

xy

forallorientationscompared tothe single-interfacemodel.For

θ

equal to0° (90°), interfacialjumps (Fig. 15d)are absent and the effectiveYoung’s modulus andPoisson’s ratiovaluesare equal in bothmodels.

Degradationis attributedto the crack-like features of the in-terfaceandisnot expectedwhenasingle-interfacemodelisused orunder a perfect interface condition. In thesesituations, inclu-sion neutrality is expectedat the previously defined neutral ori-entations.Moreover,theminimumanglesforthedouble-interface model, indicated by red marks in Fig. 15a, and the degradation rangearenotonlyafunction oftheinterface stiffnessvalues,but arealsoaffectedbygeometricalproperties(e.g.,fiberlength, num-beroffibers,andtheirspatialarrangements).

5. Micromechanicalanalysis

Mesh-independent fiber models are a necessary tool for the modelingandsimulationofcompositeswithmanyfibers(Fig.16). Tothisend,a detailedanalysis,takingintoaccount thefiber vol-umefractionasavariable,ispresentedinthisconcludingsection. Estimatesobtainedfromcommonlyusedmicro-mechanicalmodels summarizedinAppendixBarecontrastedtoresultsobtainedwith thedimensionally-reducedmodel.Theoccurrenceoffiber neutral-ityisalsoinvestigated.

Weconsiderequally-shapedfiberswithlengthlf=0.2Land

di-ameter d_f=0_.007L that are n times stiffer than the matrix ma-terial(n=10, 100, 500). To comparethe numericalresultswith thoseobtainedwithanalyticalmicromechanicalmodels,fibersand matrixare perfectlybonded.The simulationboxisdiscretized us-ing a uniform grid of 31× 31× 31 trilinear hexahedral elements, andeachfiberissubdividedinto20equally-spacedsegments.The mesh-convergence study in Fig. 17 confirms that the employed discretization is sufficiently accurate. As mentioned in the intro-duction, discontinuities in the stress field across a fiber are ne-glected indimensionally-reduced models to favornumerical effi-ciency.Theerrorintroducedwiththisapproximationisnot negligi-bleforhighstiffnesscontrastvaluesnandlargefibervolume frac-tions

υ

f.Insuchcasesaslowerconvergenceisobserved(Fig.17c).

Distributions ofhomogeneously dispersed periodic fibers are ob-tained with a random sequential adsorption algorithm described inAppendixC.Inthesimulations,fibervolumefractionsupto15% (approximately22000fibers)areconsidered.Twosample distribu-tions with approximately4700 alignedand randomly distributed discretefibersareshowninFig.16yieldingcompositeswith trans-verselyisotropicandnearlyisotropicmechanicalresponses.

Fig. 15. Effective (a) Young’s modulus, (b) shear modulus and (c) Poisson’s ratios of a two-dimensional composite with imperfect interface ( K bt1 = 10 0 0 N / mm 2 , K bt2 =

500 N / mm 2 , E

f /E m = 10 and νm = 0 . 4 in plane strain). The single-interface model shows neutrality for θ= θn ( θn ≈ 50.7 ° from Fig. 7 b), while the double-interface model

predicts degraded effective mechanical properties compared with intact matrix properties. (d) Fiber interfacial slips for θ= θn using the double-interface model. (For inter-

pretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 16. Two realizations of a composite with approximately 4700 fiber (length l f = 0 . 2 L and diameter d f = 0 . 007 L ) in an L × L × L ( L = 1 mm) periodic simulation box resulting

in a fiber volume fraction υf = 3% : (a) aligned and (b) randomly distributed fibers.

Fig. 17. Effective Young’s modulus along the x -axis for a transversely isotropic composite with aligned fibers for different mesh densities ( nn indicates the number of mesh nodes) and fiber volume fractions υf . In the calculations reported in Section 5 ( d f = 0 . 007 L and l f = 0 . 2 L ), the finest discretization with 31 × 31 × 31 trilinear hexahedral

Fig. 18. Effective Young’s modulus along the x -axis for a transversely isotropic composite with aligned fibers (a-c) and a nearly isotropic composite with randomly distributed fibers (d-f) at different fiber volume fractions υf ( d f = 0 . 007 L and l f = 0 . 2 L ).

Fig. 19. Effective transverse Young’s moduli ( E c

y = E cz ) for a transversely isotropic composite with aligned fibers at different volume fractions υf ( d f = 0 . 007 L and l f = 0 . 2 L ).

Young’s moduli. Fig. 18showsthe homogenizedYoung’s modulus

Ec

x asafunctionofthefibervolumefraction

υ

ffordifferent

mate-rial propertiesanddistributionsoffibers;inthefigure,numerical resultsarecomparedwithrule-of-mixtures(RoM)andHalpin–Tsai estimates.AsshowninFigs.18a–cforatransverselyisotropic com-posite,duetotheslendernessofthefibers(aspectratiolf/df≈ 30),

numericalresultsaregenerallyinbetteragreementwithRoM pre-dictionsespecially forfibers withrelatively lowYoung’s modulus. ThelongitudinalYoung’smodulusincreaseswiththefibervolume fraction

υ

f and is hardly influenced by the matrix Poisson’s

ra-tio

ν

m. Althougha similar pattern is observed inFigs. 18d–f for

randomly distributedfibers,the resultingYoung’smoduli are sig-nificantlylower duetofiberorientationwithrespecttothe load-ingdirectionandareingoodagreementwiththesimplemodified VoigtestimateprovidedbyPan(1996).

InatransverselyisotropiccompositesimilartothatinFig.16a, Young’s moduli Ec

y and Ezc (=Eyc) are also enhanced for a

ma-trixmaterialwithnon-zeroPoisson’sratio.However, asshownin Fig.19,thevalueofthetransverseYoung’smodulusdoesnot gen-erallyagree withRoMestimatesandthisisattributedto the un-realistic uniform stress assumption and the absence of Poisson’s

contractioneffects(HullandClyne,1996).Foranearly incompress-ible matrix,RoMpredictions clearlyunderestimatethetransverse Young’smoduliofthetransverselyisotropiccomposite.

Poisson’sratios. Fig.20showstheeffectivePoisson’sratios ofthe transverselyisotropiccompositeforvariousfibervolumefractions. SimilartothetransverseYoung’smoduli,Poisson’sratio

ν

c

zy

(

=

ν

yzc)

showsanincreasing trendandaconsiderabledependenceon the matrixPoisson’s ratio.Thisdependencehowever isnot predicted by RoM estimates,which clearlyoverestimatesthe numerical re-sults. Conversely,

ν

c

yx

(

=

ν

zxc) decreases, in good agreement with

RoM,andtheresultsarenotinfluencedbythematrixPoisson’s ra-tio.Finally,asshowninFig.20c(andalreadyreportedinFig.9cfor afixedvolumefractionwith

θ

z=0◦),

ν

xyc (=

ν

xzc) remainsequalto

thematrixPoisson’s ratioirrespective offiberdensity,a property whichisnotpredictedbyRoMestimates(Eq.(B.5)).

Figs.21and22showresultsforthehomogenizedeffective Pois-son’s ratio ofthe nearly isotropic composite asa function of

ν

m

and

υ

_f.ValuesoftheeffectivePoisson’sratio

ν

c

iso=

λ

c

Fig. 20. Effective Poisson’s ratios of a transversely isotropic composite with aligned fibers for different volume fractions ( E f /E m = 100 , d f = 0 . 007 L and l f = 0 . 2 L ).

Fig. 21. Effective Poisson’s ratios νc

iso of a nearly isotropic composite with randomly

distributed fibers for different volume fractions υf (matrix Poisson’s ratios νm ,

Ef /E m = 100 , d f = 0 . 007 L and l f = 0 . 2 L ).

where

λ

c _and

_μ

c _{are the effectiveLamé parameters averaged in}

thethree directions, are computed assuming an isotropic behav-iorofthecomposite(Nemat-NasserandHori,1998).FromFig.21, it is interesting to notice that for a matrix with Poisson’s ratio

ν

m≈ 0.25, the effective Poisson’s ratio

ν

_isoc of the homogenized

composite is insensitive to the fiber volume fraction (the inclu-sions are neutral in terms of Poisson’s ratio). For

ν

m<0.25, an

increased incompressibility of the composite is observed; for a matrixwith

ν

m>0.25, anincreasedcompressibility,althoughless

pronouncedthan theincreaseinincompressibility,isobserved.In bothcases,valuesoftheeffectivePoisson’sratiotendstowardthe value 0.25by increasing the fiber volume fraction

υ

f. As shown

inFig. 22for differentfiberstiffnesses, a similar observationcan bemadeintermsofthemeananisotropicPoisson’sratios

ν

c

a

(av-eraged overanisotropic Poisson’s ratios

ν

xy,

ν

xz,

ν

yx,

ν

zx,

ν

zy and

ν

yz,fromEq. (A.2)).Duetotherandomdistributionofthefibers,

anisotropic Poisson’s ratios only slightly deviate from the mean value.Theseobservationsareinagreementwithpreviousfindings by Das andMacKintosh (2010) who used a mean-field approach forthetheoreticalinvestigationofthemicromechanicsofisotropic compositesconsistingofrandomlydistributedstiff fibers.Bytaking intoaccountthereinforcingeffectofindividual fibers,thecurrent numericalhomogenization procedure validatesaveraging theories andneutralityintermsofPoisson’sratioforacompositewith ran-domlydistributedlinearelasticfibersembeddedinamatrixwith Poisson’s ratio

ν

m=0.25 or a nearly incompressible matrixwith

ν

m≈ 0.5.

Unlike mean-field homogenization, embedded reinforcement methods are not limited to composites with perfectmatrix-fiber

Fig. 22. Effective Poisson’s ratio of a nearly isotropic composite with randomly distributed fibers at two volume fractions υf and different matrix Poisson’s ratios νm , with

interfaces.Variousdistributionoffiberswitharbitraryshapesand interface properties can be considered for the evaluation of the compositehomogenizedmechanicalproperties.

6. Concludingremarks

It iswellknownthat fiberorientation playsan importantrole in the mechanical response of composites. The common under-standing is that fibers improvethe mechanicalproperties of the matrixinwhichthey areembedded, andthatcertain fiber orien-tationsarebetterthanothersinsomemeasure.Thisisnotcorrect. The mostunfavorable fiberorientation representsa state offiber neutrality, asituationinwhichafiberdoesnotperturbthestress field andthereforehasnoinfluenceon themechanicalproperties ofthecomposite.Fiberneutralitycanbeseenasthegeneralization oftheconceptofinclusionneutralitythatwasintroducedfor zero-thickness rigid inclusions by Wang et al. (1985). Here, we have demonstrated that inclusion neutrality holds also for thin linear elastic fibers withimperfect matrix-fiber interfaces, with an im-mediategeneralizationtotheperfectinterfacecase.

Similar to rigid line inclusions, slenderness is a fundamental ingredient for the occurrence of neutrality. Neutral fibers are in fact not reported in existing studies where composites with rel-atively large fiberdiameters were studied. In thesecases the ef-fective stiffnessofthecompositeincreaseseventhoughfibersare orientedatthemostunfavorableangle—thisanglewouldresultin fiberneutralityifthefiberwerethin.

Acknowledgments

Theresearchleadingtotheseresultshasreceivedfundingfrom the European Research Council under the European Union’s Sev-enth Framework Programme (FP7/2007-2013) / ERC Grant agree-mentno.617972.

AppendixA. Computationofeffectiveproperties

Toextracttheeffectiveelasticcoefficientsofthefiber-reinforced composite,weemploythecomputational homogenizationscheme describedbyBergeretal.(2005).AssuggestedbyXiaetal.(2003), homogeneous boundary conditions lead to over-constrained pre-dictionsoftheeffectiveproperties.Instead,theperiodicboundary conditions

u k

i − uli=

ε

¯i jL j, (A.1)

prescribed at theboundary ofthe simulation boxwith

ε

¯i j

(

i,j=

x,y,z

)

thecomponentsoftheappliedstraintensor,yieldmore ac-curate estimations. Edge indexes k and l correspond to two op-posite edgesofthesimulationboxwhereperiodicstrainsare im-posed,andLj isthesize ofthe domaininthejthdirection. Inall

thenumericalresultspresentedinthispaper,aunitcube simula-tionboxisadopted(Fig.3).

Although21independentconstantsare definedforthegeneral anisotropiccase(Nemat-NasserandHori,1998),hereonlythe spe-cificengineeringconstantsineach directionareofinterest.Three Young’s moduli (Ec

x, Ecy, Ezc) relating tensile stresses and strains,

three shearmoduli (Gc

xy, Gcyz, Gcxz) relatingshearing stresses and

strains, andsix Poisson’sratios

ν

c

i j (

ν

xyc,

ν

xzc,

ν

yxc,

ν

yzc,

ν

zxc,

ν

zyc)

representingcontractioninthej-directionafterthe applicationof atensileloadinthei-directioncan beextracted.Thesequantities may be expressed in termsof the coefficientsof the compliance constantsmatrixas E c x = 1 C 11, E c y= 1 C 22, E c z= 1 C 33, G c xy= 1 C 44 , G c yz= 1 C 55 , G c xz= 1 C 66 ,

ν

c xy=− C 21 C 11,

ν

c xz=− C 31 C 11,

ν

c yx=− C 12 C 22 ,

ν

c yz=− C 32 C 22 ,

ν

c zx=− C 13 C 33,

ν

c zy=− C 23 C 33. (A.2) The compliance constants matrixC is evaluated by inverting the elasticcoefficientmatrixDasdescribedinMalagù etal.(2017).

AppendixB. Analyticalmicromechanicalmodels

Forthespecialcaseofatransversely isotropiccompositewith alignedfibersofthetypeshowninFig.16a,simpleanalytical mi-cromechanicalmodelsareavailable.

Thewellknownruleofmixturesiswidelyusedto predictthe mechanical propertiesof elastic composites. The composite stiff-ness is approximated as a weightedmean of the moduli oftwo components(HullandClyne,1996).Theruleofmixturesdoesnot consider geometrical details such asthe fiber aspect ratio. With referenceto thecaseof a compositewith fibersthat are aligned withtheglobalx-axis,the effectivelongitudinalYoung’smodulus

E c

x=

(

1−

υ

f

)

E m+

υ

fE f (B.1)

isexpressed asafunction ofthe fibervolume fraction

υ

_f, where itis assumedthat strainsin thedirectionof thefibers are equal inthematrixandthefiber(this assumptionisvalidforperfectly bondedfibers withhighaspect ratio). Eq. (B.1) isreferred to as theupperboundmodulus(Voigtmodel).

Usingthe inverserule ofmixtures andassuming that stresses inthedirectionnormal tothe fibersare equalinthe matrixand thefiber,thetransverseYoung’smoduli

E c y=E zc=



_υ

f E f + 1−

υ

f E m



−1 . (B.2)

This expression is known as the lower-bound modulus (Reuss model).

Accordingtothesemi-empiricalestimateoftheYoung’s modu-lusproposedbyHalpin(1969),

E c x = E m 1+

ξηυ

f 1−

ηυ

f (B.3) with

η

= E f/E m − 1 E f/E m +

ξ

and

ξ

= 2l f/d f. (B.4)

Thisexpressionissuitable forshortfibersandtakesintoaccount thenon-uniformdistributionofstrainsandstressesinthe compos-ite.PredictionsofEq. (B.3)lie withintherangeofthelowerand upperboundmoduli.

Similar formulas can be derived for the Poisson’s ra-tios(HullandClyne,1996):

ν

c xy=

(

1−

υ

f

)

ν

m+

υ

f

ν

f,

ν

c yx=[

(

1−

υ

f

)

ν

m+

υ

f

ν

f] E y E x, and

ν

c yz=1−

ν

yx− E y 3K with K =



_υ

_f K f + 1−

υ

f K m



−1 , (B.5) whereKm andKf arethebulkmodulioftheconstituents.

Foracompositewithrandomlydistributedfiberswithisotropic behavior,thesimpleruleofmixturesestimateyields

E c x,y,z=



1−

υ

f 2

π



E m+

υ

f 2

π

E f, (B.6)

whereaprobabilitydensityfunctionisusedtospecifytherandom fiberorientationasproposedbyPan(1996).