Effective behaviour of porous ductile solids with a non-quadratic isotropic matrix yield surface

ContentslistsavailableatScienceDirect

Journal of the Mechanics and Physics of Solids

journalhomepage:www.elsevier.com/locate/jmps

Effective behaviour of porous ductile solids with a non-quadratic isotropic matrix yield surface

Lars Edvard Blystad Dæhli

^a,∗

, Odd Sture Hopperstad

^a

, Ahmed Benallal

^b

aStructural Impact Laboratory (SIMLab), Department of Structural Engineering, Norwegian University of Science and Technology (NTNU), NO-7491 Trondheim, Norway

bLMT, ENS Paris-Saclay/CNRS/Université Paris-Saclay, 61 Avenue du Président Wilson, 94235 Cachan, France

a rt i c l e i n f o

Article history:

Received 21 November 2018 Revised 16 April 2019 Accepted 20 May 2019 Available online 21 May 2019 Keywords:

Numerical limit analysis Third deviatoric stress invariant Unit cell modelling

Porous plasticity

a b s t ra c t

Inthisstudy,weexaminethemacroscopicyieldingofisotropicporousductilesolidshavinga matrixyieldfunctiondependentonthesecondandthirddeviatoricstressinvariants.Numer- icallimitanalysesusingathree-dimensionalfiniteelementmodelofahollowspherewitha Hershey-Hosfordmatrixyieldfunctionareconductedfordifferentshapesofthematrixyield surfaceandporositylevels.Thesenumericalresultsarethenusedtoelucidatefirst-orderef- fects ofthethirddeviatoric stressinvariant onthemacroscopicyielding andfurtherused asreferencedatatoassesstheperformanceof twoporous plasticitymodelsthatincorpo- rateeffectsofthethirddeviatoricstressinvariantusingtheisotropicnon-quadraticHershey- Hosfordyieldfunction.Thefirstmodelisderivedfromanupper-boundlimitanalysisofthe hollowsphererepresentativevolumeelementusingtheGurson-Ricetrialvelocityfield,but witharathergeneralisotropicmatrixyieldfunction.Thesecondmodelisasimple,heuristic extensionoftheGursonmodelincorporatingtheequivalentstressmeasureoftheHershey- Hosfordyieldfunction.

Fromthenumericallimitanalyses,itisfoundthatthecontoursofthemacroscopicyield surfaceinthedeviatoricplanetransformfromthehexagonalshapeoftheunderlyingmatrix yieldsurfacetoaroundedtriangularshapethatconvergestothecircularshapeoftheGur- sonmodelasthemacroscopicstresstriaxialityratioincreases.Thisshapetransformationis dependentupontheporositylevel.Theupper-boundmodel wasfoundtobeinverygood agreementwiththenumericaldataforallstressstates,shapesofthematrixyieldsurface, andporositylevels.Theheuristicmodelprovidesgoodpredictionsforlowandmoderatelev- elsofporositypertinenttoductilefracture,butthepredictionsdeterioratewhenthestress triaxialityratioandtheporositylevelincrease.

Wealsoaddresstheissueofhowrepresentativethesphericalunitcellisforthedescrip- tionofrealporoussolids.Tothatend,wemakecomparisonsbetweenaspace-fillingrepre- sentativevolumeelementintheformofacubicunitcellwithacentricsphericalvoidand thehollowspheremodel.Theseresultsshowthatthehollowsphere modelgenerallypro- videsslightlyhighervaluesfortheyieldlimits.Theshapeoftheyieldlociissimilarforthe twomodelsinthecaseofnon-quadraticmatrixyieldsurfaces,whilethecubicmodelgives adifferentshapeoftheyieldlociforlowandintermediatestresstriaxialityratioswhenthe vonMisesyieldfunctionisusedforthematrix.

© 2019TheAuthors.PublishedbyElsevierLtd.

ThisisanopenaccessarticleundertheCCBY-NC-NDlicense.

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

∗ Corresponding author.

E-mail address: lars.e.dahli@ntnu.no (L.E.B. Dæhli).

https://doi.org/10.1016/j.jmps.2019.05.014

0022-5096/© 2019 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license.

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

1. Introduction

Theyieldingandplasticflowofmanymetalalloysaredependent onthethirddeviatoricstressinvariant (J₃).Theyield surfacehasregions ofhigher curvatureandis thus anon-quadratic function ofthe stress components.The shapeof the yield surfacearisesfromthe physical mechanismsof plasticslip ina polycrystal material anddependsmarkedly on the crystallographictextureofthematerial.However,alsomaterialswithrandomtextureexhibitanon-quadraticyieldsurface.

Forinstance,FCCmaterialssuchasaluminiumalloyshavebeenshowntobemoreaccuratelydescribedbyanon-quadratic yieldsurfaceeveninthecasewhenthematerialisisotropic(Hosford,1972,1996;LianandChen,1991).

Porousplasticitymodelsaretraditionallydevelopedbasedonarepresentativevolumeelement(RVE)consistingofmatrix materialandasphericalvoid,wherethematrixmaterialisgovernedbythevonMisesyieldfunction.AsthevonMisesyield functiondependsonlyonthesecond deviatoricstress invariant(J₂),thereisnomicroscopicdependenceofJ₃ throughthe matrixconstitutive formulation.However,aslightJ₃ dependencestill appearsonthemacroscopiclevel,i.e.inthehomog- enizedresponseoftheRVE(Benallaletal.,2014;CazacuandRevil-Baudard,2015;Cazacu etal., 2013;LeblondandMorin, 2014). Thisstems fromthe stress heterogeneityimposed by the voidedstructure ofthe RVE.Thismacroscopic J₃ depen- dencewasindeedeliminatedintheclassicderivationoftheporousplasticitymodelbyGurson(1977)duetotheexclusion ofhigher-orderterms inthe Taylorseriesexpansion ofthe microscopic plasticdissipation function (Benallaletal., 2014;

Cazacuetal.,2013;LeblondandMorin,2014).

3DyieldsurfacesforaporousmaterialwithvonMisesandTrescamatrixbehaviourwerepresentedinRevil-Baudardand Cazacu (2014a). The yieldsurfaces were shownto be affected by J₃ forboth matrixdescriptions, howevermost severely forthe Trescamatrix behaviour,andit wasfound that the yieldsurfaces were centro-symmetric.Also, theporosity was seento be ofkey importance forhow quicklytheyield locichanged shape withincreasing levels ofhydrostatictension.

Thisstudydemonstrates theintricatecouplingbetweenthestressinvariants, thematrixbehaviour andtheporosity level.

Isotropicyieldcriteriaformatrixmaterialsincorporatingfirst-ordereffectsofJ3 basedonlimitanalysisusingtheGurson- Ricetrialvelocityfield(Gurson,1977;RiceandTracey,1969)havebeenconsideredinrecentstudiesbyCazacuetal.(2014b), Soare(2016)andBenallal(2017,2018).Cazacuetal.(2014b)adopttheTrescayieldcriterionforthematrixmaterialandde- riveananalyticalcriterionforaporoussolidunderaxisymmetricloadings.Soare(2016)employsthenon-quadraticisotropic yield criterion first proposed by Hershey (1954) and later by Hosford (1972), which we henceforth will refer to as the Hershey-Hosfordyieldcriterion. Benallal(2017, 2018)considers ageneralisotropic yieldcriterion forthematrix material, whichincludestheHershey-Hosfordyieldcriterion.ThestudiesbySoare(2016)andBenallal(2017,2018)arenot restricted toaxisymmetricloadingconditionsandresultinayieldsurfacefortheporoussolidthatdramaticallychangesshapeinthe ^{-plane asthe stress statechanges from}predominantly deviatoricto predominantlyhydrostatic. More specifically, when thestress stateispredominantly deviatoric,the yieldsurfaceoftheporoussolid resemblesthat oftheunderlyingmatrix material.Whenthemagnitudeofthehydrostaticstressincreases,thecontoursoftheyieldsurfaceoftheporoussolidfor givenhydrostaticstressstatesdisplayaroundedtriangularshapeinthedeviatoricplane– aneffectthatisexplicitlyshown inthepapersbySoare(2016)andBenallal(2017).Moreover,theyieldsurfacedisplayscentro-symmetry,whichisageneral propertyforaporousplasticmaterialwherethematrixyieldfunctionisanevenfunctionofthestressstate(Cazacuetal., 2019).

Apart from the analytical work utilizing limit analysis, there are also a number of studies that have been devoted to numerical limit analyses. Trillat and Pastor (2005) made a thorough comparison between static (lower bound) and kinematic (upper bound) numerical limit analyses for a hollow sphere RVE and the Gurson model. Theyconclude that Gurson’s criterion is indeeda good approximation. FE limit analyses have also recently been performed by Madou and Leblond (2012) where both ellipsoidal (prolate and oblate) and cylindrical geometries were considered and used to assess their analytical porous plasticity model. The effects associated with an assemblage of voids were studied by Fritzenetal.(2012).Inthatstudy,FEsimulationsofporousmaterials comprisingarandomspatialdistributionofspherical voidsinanelasto-perfectplasticmaterialgovernedbyJ₂ flowtheorywereperformed.Theauthorsconsideredvariousspa- tialconfigurationsandevaluatedthestatisticalpropertiesoftheunitcellsimulations.Interestingly,theirresultsshowthat anapproximationbasedonasinglecenteredvoid,asisthecaseforthehollowspheremodel,givesfairlygoodagreement toarandomdistributionofvoidswhentheporosityisratherlow,whichintheirstudycorrespondstoaporosityof1%.This isaratherrealistic volumefractionofpotentiallyvoid-nucleatingparticles formanystructuralmetals.Moreover,their re- sultsshowthatthepredictionsoftheGursonmodelarelessaccurateforhighporositylevels(upto f=0.3intheirstudy).

ThisiscorroboratedbytheresultsinCazacuetal.(2014b)foraTrescamatrixusingporositylevelsbetween f=0.001and f=0.04.Alsotheparticleshapecanaffectthenumericalresults,whichwasstudiedrecentlybyKeralavarma(2017)using bothcubicandsphericalvoids innumericalFElimitanalyses.The studyshowsthattheRVEwithacubicvoidlowersthe resultingplasticlimitload;however,onlyslightly,andtheresultswerequalitativelysimilar

Othertypesofmatrixyieldfunctionshavealsobeenexaminedintermsofnumericallimitanalyses.Guoetal.(2008)em- ployed a Drucker-Prager matrix yield function and compared the FE results to their analytical criterion. In a study by Pastoret al.(2010), numericallimit analysisofthe hollow sphere modelwith CoulombandDrucker-Prager matrix yield surfaceswasundertakenandcomparedtoexistinganalyticalmodels.Thoré etal.(2011)usedbothvonMisesandDrucker- Pragermatrixformulationsandtheirstudyconfirms,fromanumericalpoint-of-view,theJ₃dependenceofthemacroscopic yieldsurfacefortheporousplasticsolid whenthematrixisgovernedbythevon Misesyieldcriterion.Similarresultsfor avonMisesmatrixwereobtainedbyAlves etal.(2014)withanFEmodelcorresponding toacubicRVEwithacentered,

spherical void. ATresca matrix formulationwas employed in axisymmetric FE calculationsby Cazacu etal. (2014b) and further extended tofull three-dimensional FEsimulations in Cazacu etal. (2014a) utilizinga cubic RVEwitha centered, sphericalvoid.Pastoretal.(2015)performednumericallimitanalyses(bothstaticandkinematic)ofaporousmaterialcom- prisingellipticcylindricalvoidsembeddedinaCoulombmatrixmaterial.Tothebestknowledgeoftheauthors,nostudies havesofarconsideredFElimit analysisusingtheHershey-Hosfordyieldcriterion.Thisyieldcriterionisofgreatrelevance forpolycrystallinesolidsthatarefrequentlyusedforindustrialapplications,e.g. foraluminiumalloysthat areincreasingly employedbytheautomotiveindustryandinprotectivestructures.Thisservesasakeymotivationforthecurrentstudy.

ItwasshownbyBarlat (1987)howtheshapeoftheyieldsurfaceforisotropicsolidsaffectstheformabilityofasheet materialdeformedunderplanestressconditions.Marciniak-Kuczy´nskianalyses (MarciniakandKuczy´nski,1967)were con- ducted withtheHershey-Hosford yieldcriterion and itwasfound that the formabilitydecreased dramaticallyforbiaxial stressstateswithincreasingcurvatureoftheyieldsurface.Recently,theinfluenceoftheyieldsurfaceshapeonductilefail- ureby strainlocalizationwasinvestigatedbyDæhlietal.(2017b).Tothisend, numericalsimulationswithaunit cellthat resemblesan imperfectionband wereperformed, adoptingtheHershey-Hosfordyieldcriterion forthematrixmaterial.In addition,aheuristicextensionoftheGursonmodelthatincorporatesJ₃dependenceviatheHershey-Hosfordyieldfunction wasproposed andevaluatedagainstunitcellsimulations,andthenusedinlocalizationanalyses accordingtothetheoreti- cal frameworkestablishedbyRice(1976).Itwasdemonstratedbybothmethodsthat theshapeoftheyieldsurfaceofthe underlyingmatrixmaterialhasapronouncedinfluenceontheonsetofstrain localization,andthetrendisthat thestrain to failure decreases whenthe curvatureof theyield surfacebecomeshigher. Thishighlights theimportance ofincluding first-ordereffectsofJ₃ inporousplasticitymodels andsubstantiatestherelevanceofmakingaproperassessmentofsuch models.

Inthiswork,finiteelementsimulationsofanRVEintheformofahollowsphereareconductedtoevaluatethemacro- scopicyieldsurfaceofaporoussolid.Hence,plasticyieldingisexpressedintermsofstresscomponentsthatareaveraged overthevolumeoftheRVE.PerfectplasticmaterialbehaviourgovernedbytheHershey-Hosfordyieldcriterionandtheas- sociativeflowruleisadoptedforthematrixmaterialtoevaluatethemacroscopicyieldsurfacenumerically.Theexponentof theHershey-Hosfordyieldfunctionisvariedtoinvestigatetheinfluenceofthematrixyieldsurfacecurvatureontheshape ofthemacroscopicyieldsurface.Itisfoundthattheshapeofthemacroscopicyieldlocusinthe^{-planedependsmarkedly} onthemacroscopicstresstriaxialitywhentheunderlyingmatrixmaterialisJ₃dependent.Thenumericalresultsarefurther usedtoassesstheperformance ofthetwoporousplasticitymodelsbasedontheHershey-Hosfordyieldfunctionrecently proposedbyBenallal(2017)andDæhlietal.(2017b).ItisfoundthattheporousplasticitymodelbyBenallal(2017),which isderived froma rigorousupper-bound limitanalysis, isin verygoodagreement withtheRVEsimulations.The heuristic porousplasticitymodelusedbyDæhlietal.(2017b)isfoundtogiveareasonablerepresentationoftheyieldsurfaceforthe porositylevelstypicallyfoundinmetalalloys,butitislessaccuratethantheupper-boundmodelforhighstress triaxiality ratiosandhighporositylevels.

ThenotationusedthroughoutthepaperisgiveninSection2.Theconstitutivemodelofthematrixmaterialispresented inSection3,whileSection4brieflypresentsthetwoporousplasticitymodels.Section5dealswiththeformulationofthe finite elementmodel fortheRVE.The resultsare presentedanddiscussed inSection 6.Concluding remarks aregivenin Section7.

2. Notation

We adoptaCartesian coordinatesystem(x₁,x₂,x₃) withorthonormalbasevectors(e₁,e₂,e₃), andexpressall tensor componentsinthiscoordinatesystemifnothingelseisstated.TheEinsteinsummationconventionisusedthroughoutthe paper,unlessspecificallystatedotherwise.

Themicroscopic(orlocal)stressandstrainratetensorsofthematrixmaterialaredenoted

σ

=

σ

i je_ie_jand

ε

˙=

ε

˙i je_i e_j,respectively.The invariantsofthemicroscopic stresstensor usedinthiswork arethehydrostaticstress (

σ

h), thevon Misesequivalentstress(

σ

eq^vm),andthedeviatoricangle(

ω

^{).Theseinvariantsareexpressedas}

σ

h=

σ

kk

3 (1)

σ

eq^vm=

3

2

σ

i j

σ

i j (2)

ω

=1 3arccos

27 2

det

( σ

) ( σ

eq^vm

)

³

(3) where

σ

=

σ

i j−

σ

h

δ

i j

e_ie_jisthemicroscopicdeviatoricstresstensorand

δ

ij istheKroneckerdelta.

Themacroscopicstressandstrainratetensorsaredefinedbyaveragingtheirmicroscopiccounterpartsoverthevolume VoftheRVE,i.e.

=

i je_ie_j= 1 V

V

σ

^{dV (4)}

E˙ =E˙_{i j}e_ie_j= 1 V

V

ε

˙ ^{dV (5)}

whereVcomprisesthevolumeofthevoidVvandthevolumeofthematrixVm,viz.

V=Vv+Vm (6)

Theinvariantsofthemacroscopicstresstensoraregivenby

h=

kk

3 (7)

^vmeq =

3

2

_{i j}

_{i j} ⁽⁸⁾

= 1 3arccos

27 2

det

(

) (

eq^vm

)

³

(9)

wheretheprimeagainreferstothedeviatoricpartofthetensor.Wenotethatthemacroscopicdeviatoricangle(^)specifies theanglebetweentheprojected1-axis(m₁)andthestresspoint(P)inthedeviatoricplane,asillustratedinFig.1(b).The definitionof

ω

^forthemicroscopicstresstensoriscompletelyanalogous.

Wewillalsoreferextensivelytothemacroscopicstresstriaxialityratio,whichisdefinedas T=

h

^vmeq

(10)

The finite element (FE) limit analyses will be conducted under proportional macroscopic loading. To this end, the macroscopic stress triaxiality (T) and the macroscopic deviatoric angle (^{) are kept constant throughout the loading} sequence.

Thestresstensordecomposesintoadeviatoricandahydrostaticpart.Themacroscopicprincipalstresscomponentsmay consequentlybewrittenonageneralformusingthedimensionlessparametersTand^,viz.

1=

₂

3cos

( )

^+T

eq^vm (11a)

2=

₂

3cos

−2

π

3 +T

eq^vm (11b)

3=

₂

3cos

⁺²

π

3 +T

eq^vm (11c)

Wenote thatthesearegenerallynottheorderedprincipalstresscomponents,sincetheirorderingisexclusivelyrelated tothedeviatoricangle^{.However,inthecurrentcontext,weemployisotropyandthedeviatoricanglecanthenbere-} strictedto0^◦≤≤60^◦.Thus,theorderingofthemacroscopicprincipalstressescorrespondsto1≥2≥3inthecurrent study.Fig.1providesagraphicalinterpretationofthestressstateparametersusedinthesequel.

Fig. 1. Illustration of an arbitrary stress point ( P ) depicted in ( a ) the principal stress space and ( b ) the deviatoric plane. The unit vector e his directed along the hydrostatic axis and the three vectors m iare the projected unit base vectors in the deviatoric plane.

3. Matrixconstitutivemodel

TheelasticresponseofthematrixmaterialisgovernedbythehypoelasticformulationofthegeneralizedHooke’slaw,for whichtheelasticmodulus(E)andthePoissonratio(

ν

^{)arelistedinTable1.WesetthePoissonratiototheratherhighvalue}

of

ν

=0.49inthe numericalcalculationsto mimicaclose-to-incompressible materialresponse alsointheelasticdomain withoutintroducingnumericalstabilityproblems(MadouandLeblond,2012;2013;Teko˘gluetal.,2012).Thisensuresthat an incompressiblevelocity fielddevelops soonerin theFEcalculations, whichis aprerequisite fortheupper-bound limit analysisperformedbyGurson(1977).

TheplasticresponseofthematrixmaterialisgovernedbytheHershey-Hosfordnon-quadraticisotropicyieldfunction

φ ( σ )

=

σ

^eq

( σ )

−

σ

0≤0 (12)

wheretheequivalentstressisgivenby

σ

^eq=

₁

2

( σ

I−

σ

II

)

^m+

( σ

I−

σ

III

)

^m+

( σ

II−

σ

III

)

^m

m¹

=

σ

eq^vmg

( ω )

⁽¹³⁾

The orderedprincipalstress componentsare denoted

σ

I≥

σ

II≥

σ

III andm representstheyieldsurfaceexponentthat dic- tatestheshapeoftheyieldsurface. Theequivalentstressmaybe writtenasa productofthevonMisesequivalentstress and a function g(

ω

^{) that depends}non-linearly on the deviatoricangle, asindicated by the last equalityin Eq. (13).We note that exponentvalues within theranges 1≤m≤2 and4≤m<∞give yield surfaceslying in-between the von Mises (m=2,4)andtheTresca(m=1,∞)yield surfaces.However, thehigherrangeofexponentsismoreconvenientandmost frequently used. The exponent values m=6 and m=8 are often associated with BCC and FCC materials, respectively (Hosford,1996).

Theyieldstressinuniaxialtensionisdenoted

σ

0andisconstantinthisstudytoconformwiththeframeworkofclassic limitanalysis,whichhingesonperfectplasticmaterialbehaviour.Further,weadopttheassociatedflowruleandtheplastic strainratetensorisaccordinglygivenby

ε

˙^p=

λ∂φ

^˙

⁽ σ )

∂ σ

⁽¹⁴⁾

where

λ

˙ ≥0isthe non-negativeplasticloading parameter.The matrixyield stressandthe yieldsurfaceexponentvalues employed in thisstudy are listed in Table1. Fig. 2shows the matrixyield surfaces corresponding to the differentyield surfaceexponentsusedinthenumericalcalculations.

AUMATsubroutineavailablethrough themateriallibraryattheSIMLab (StructuralImpactLaboratory)groupatNTNU (NorwegianUniversityofScienceandTechnology)wasusedinthenumericalsimulations.Theplasticcorrectionsoftheelas- tictrialstatewere governedbyafully-implicit Backward-Eulerstress-updatealgorithm.We alsoemployedasub-stepping proceduretoensurethatsufficientlysmallstrainincrementswereusedintheplasticiterations.Thestrainincrementssent

Table 1

Overview of the constitutive parameters used for the matrix material in the numerical analyses.

Material parameter E [MPa] ν σ0[MPa] m χ Parameter value 100 0.49 1 2, 8, 20 0.1

Fig. 2. Matrix yield surfaces for the exponent values used in this study depicted in the deviatoric plane.

tothe stress-updatealgorithm were then governed by √

ε

^:

ε

≤

χσ

0/E,where

ε

^{denotes thetotal strain increment}

and

χ

^{isa sub-step thresholdparameter.In thiswork,we employ asub-step thresholdparameter value of}

χ

=0.1.This warrantsarathersmallstrainincrementtobeusedintheplasticcorrectioniterations.Wenotethattheratiobetweenthe yieldstressandtheelasticmoduluswas

σ

0/E=0.01inalltheanalysesconductedherein,asinferredfromtheparameters listedinTable1,whichcorresponds toan unusuallyhighyieldstrainfortypicalmetals.However,thisvaluehasnosignif- icanceforthecurrentstudysinceweareonlyinterestedindeterminingthelimitloadoftheunitcell.Asthelimitloadis attained,i.e.whenasufficientportionoftheunitcellbehavesperfectlyplastic,themacroscopicresponseisunaffectedby elasticitysincethestresseshavesaturated.

4. Porousplasticitymodelswithnon-quadraticmatrixyieldcriterion 4.1. Rigorousmodeldevelopedfromupper-boundanalysis

We will only make a brief presentation of the porous plasticity model in this section and the reader is referred to theoriginal article formoredetails (Benallal,2017). The model hasbeendeveloped along the samelines asthe original Gursonmodel,with an incompressible, isotropicandrigid-plastic matrixdescription,using the Gurson-Ricetrialvelocity field(Gurson,1977;RiceandTracey,1969).However,incontrasttotheoriginalGursonmodel,plasticyieldingofthematrix isnowtakentobeJ₃dependentandgovernedbythegeneralrelation

φ ( σ )

=

σ

^eq−

σ

0=

σ

eq^vmg

( ω )

−

σ

0≤0 (15) wherethe microscopic deviatoricangle (

ω

^{(x)) varieswiththe position (x) throughoutthe matrixmaterial andgenerally}

deviatesfromthemacroscopiccounterpart(^).

Themodel canformally be writtenin thesame globalexpressionasthe Gursonmodel,butit involvesfourfunctions a,R, P andQ that all dependon the macroscopicstress triaxiality T,the macroscopicdeviatoric angle^{and theporos-} ityf. These functions are generally not available in closed form, butthey can be given explicit forms by numerical in- tegration andinterpolation of the resultingdata. The effective yield surfacemay then be written on the form (Benallal, 2017)

=

a

eq^vm+P R

σ

0

2

+2fcosh

h−Q R

σ

0

−

(

¹+f²

)

=0 (16)

Inthecurrentwork,we havenotfocused onevaluatingthe functionsa, R,PandQ.Wehaveratherobtainedyield stress valuesbyimposingamultitudeofvaluesforthetriaxialityanddeviatoricangleofthestrainratetensor,whichthengives acloudofyieldstresspointsintermsofTand^.Afterwards,thecontourplotfeatureoftheMatplotliblibraryinPython wasusedtoplottheyieldsurfacebyinterpolationoftheyieldstressvaluesintheT−^space.

ItwasshownbyBenallal(2017)thatwhenT→∞,theeffectiveyieldpointtendstowardsthehydrostaticpoint

h

σ

^{0 =−}

2 3lnf

g

(

^π3

)

⁽¹⁷⁾

AsimilarresultisobtainedwhenT→−∞,whichleadsto

h

σ

^{0 =}

2 3lnf

g

(

⁰

)

⁽¹⁸⁾

WenotethatsincetheHershey-Hosfordyieldfunctionissymmetricwithrespectto

ω

=^π₆,wehavethatg(⁰)=g(^π3)^and consequentlythat themagnitudeofthehydrostaticstressis equalforT→± ∞.Forpure shearloadingwithh=0,itis foundthattheequationoftheyieldlocusinthe^{-planeisgivenby}

^vmeqg

()

⁼

(

^1−f

) σ

0 (19)

whichisequaltothematrixyieldfunction (seeEq.(15))up tothesizereductionfactor1−f.However,note thatinthis caseitisthemacroscopicdeviatoricangle(^{)thatentersthenon-linearfunctiong.}

Inthiswork,wewillfocussolelyonmatrixconstitutivebehaviourgovernedbytheHershey-Hosfordnon-quadraticyield function givenin Eq.(12). However, theyield function for theporous plastic materialshown inEq. (16)applies to gen- eraltypes of isotropic J₃-dependent matrix yield functions. Henceforth, we will refer to thismodel as the upper-bound modelor the rigorousmodel. We note that rigoroushere means that the upperbound isobtained in a mathematically rigorous way using the Gurson-Rice trial velocity field, without introducing any simplifying assumptions regarding the microscopic dissipation function. We also note that the modification introduced by Tvergaard (1981) can be employed in the porous plasticity model governed by Eq. (16), which can serve as a means to enhance the predictions in some cases.

4.2. HeuristicextensionoftheGursonmodel

Theporousplasticitymodeldevelopedintheprevioussection correspondstotheexactsolutionofthekinematiclimit analysisusingtheGurson-Ricetrialvelocityfield. Anotherwayto includea J₃ dependenceintheGurson modelwaspre- sentedin Dæhlietal.(2017b)basedon thework ofDoege andSeibert(1995),wherethey propose aheuristic extension oftheGursonmodeltoincludematrixconstitutivebehaviourgovernedbytheHill48anisotropicyieldfunction.Themodel usedbyDæhlietal.(2017b)preservesthestructureoftheGursonyieldcriterion,viz.

=

^eq

σ

0

2

+2q₁fcosh

3 2q₂

h

σ

0

−

(

¹+q₃f²

)

=0 (20)

wherethemodificationsintroducedbyTvergaard(1981)havebeenincludedsincethisformismostfrequentlyadoptedin theliterature. Theonly differencebetweentheGurson-Tvergaardmodel andthisheuristic extension istheexpression for themacroscopicequivalentstresseq,whichisnowdefinedas

eq=

₁

2

(

I−

II

)

^m+

(

I−

III

)

^m+

(

II−

III

)

^m

m¹

=

^vmeqg

()

⁽²¹⁾

whereI≥II≥IIIaretheorderedmacroscopicprincipalstresses.Ingeneral,gisanon-linearfunctionofthemacroscopic deviatoricangle(^{),butitisconstantandequaltounityform}=2,4forwhichEq.(20)reducestotheGurson-Tvergaard model.Throughoutthiswork,wewillretainthestructureoftheoriginalGursonmodelandthussetq₁=q₂=q₃=1.We notethatEq.(20)isthenexactlytheGursonmodel,howeverwiththevonMisesequivalentstressreplacedbytheequivalent stressinEq.(21).

4.3. Comparisonbetweenthetwoporousplasticitymodels

TheyieldsurfacesobtainedbythetwoporousplasticitymodelspresentedinSections4.1and4.2throughEqs.(16)and (20) are compared in Fig. 3 for yield surface exponents m=2,8,20 and porosity levels f=0.001,0.01,0.1. These fig- uresshow yield lociforcontours of constanthydrostatic stress plottedinthe plane spannedby the axeseq^vmcos ^{and vm}eq sin^.

ItisevidentfromFig.3thatthemacroscopicyieldsurfacesofthetwoporousplasticitymodelsarenearlyidenticalwhen thevonMisesyieldfunctionisusedforthematrixmaterial(m=2).Thesmalldifferencesstemfromthetruncationofthe seriesapproximationofthematrixdissipationfunctionintheoriginalGursonmodelandarepredominantlyrelatedtothe macroscopicstresstriaxiality.IthasbeenshownbyLeblondandMorin(2014)thatbycalculatingthesecond-orderapprox- imationofthemacroscopicdissipation,theresultingyieldsurfaceoftheporoussolidliesinteriortotheGursonmodel.The porousplasticitymodelpresentedinSection4.1correspondstoanexactintegrationofthemacroscopicdissipationfunction, andthusincludesthiseffectofthemacroscopicstresstriaxialitythroughtheparametersa,R,PandQ(Benallaletal.,2014).

We note that also some effectsare associated withJ₃,butthey are too smallto be detected on theyield locishown in Fig.3(a)–3(c).

Whentheyieldfunctionoftheunderlyingmatrixmaterialisnon-quadratic(m=8andm=20),theshapeoftheyield surfaceofthetwoporousplasticitymodelsismarkedlydifferentfornon-zerohydrostaticstress.Therigorousupper-bound modelgivesayieldsurfacethatevolvesfromthe

π

^{/6-typesymmetry(hexagonalshape)tothe}

π

^{/3-typesymmetry(rounded}

triangularshape)asthehydrostaticstressincreases,whereastheheuristicmodelretainsthe

π

^{/6-typesymmetry(hexagonal}

shape)oftheunderlyingmatrixyieldfunctionforalllevelsofhydrostaticstress.Thisisreadilyseenfromtheyieldfunction expressedinEq.(20);whenthehydrostaticstressisconstant,theyieldfunctionisgovernedbytheequivalentstress,which corresponds tothe Hershey-Hosfordyield criterion.We note thatthe roundedtriangularshapedisplayedby the rigorous model leadsto a higheryield stress in generalizedcompression than in generalized tension forthe positive hydrostatic stresslevelsshowninFig.3.Thiseffectisoppositefornegativehydrostaticstresslevels,suchthattheyieldstressingener- alizedtensionishigherthaningeneralizedcompression,becausetheyieldfunctionsplottedinFig.3arecentro-symmetric (Cazacuetal.,2019).

The yield lociof the rigorousupper-bound modelusually lie interiorto the yield loci given by the heuristic model, andthusrepresentsatighterupper-boundsolutiontothelimit analysisinmostcases.However, exceptionstothispertain to stress states close to generalized shear (=30^◦) when the hydrostatic stress level is rather low, but different from zero. With a keen eye, this can be observed in Fig. 3(g)–3(i) for the contours corresponding to h=3.5, 2.0, and 0.9, respectively.

Theheuristicmodelandtherigorousupper-boundmodelcoincideinthe^-plane(h=0),andalsoatthehydrostatic limits(bothcompressiveandtensile)inthecaseoftheHershey-Hosfordmodelconsideredhere.However, notethatwhile theformerholdsinthegeneralcase,thelatterdoesnot,sincetheheuristicmodelalways leadstotheGursonhydrostatic limits,whiletherigorousmodelgivesthelimitsshowninEqs.(17)and(18).Theselimitsareexactforthehollowsphere RVE withgeneral non-quadraticisotropic matrix yield functions(see Benallal, 2018). In the caseof theHershey-Hosford model,wehavethatg(⁰)=g(

π

/3)=1,andthetwolimitscoincide.

Fig. 3. Comparison between yield surfaces resulting from the rigorous model (solid lines) and the heuristic model (dashed lines) for different yield surface exponents; ( a )–( c ) m = 2 , ( d )–( f ) m = 8 , and ( g )-( i ) m = 20 . The contours represent yield loci for constant levels of hydrostatic stress ( h= cnst ).

5. Finiteelementmodel

5.1. Overviewofthefiniteelementmodel

The RVEused in theFE analyses corresponds to a hollow sphere, inaccordance withthe Gurson modelformulation.

Weemploysymmetryconditionsinthethreemutuallyorthogonaldirectionsx₁,x₂,andx₃.Hence,onlyone-eighthofthe RVEis explicitlyconsidered in theFE model,asillustrated in Fig. 4. One ofthe three symmetry planes(with normalin thex₃-direction)ishighlightedwithabluecolourinFig.4,whilethetwoothersymmetryplanes(withnormalsalongthe x₁-directionandx₂-direction)arenot shown.The outersurfaceofthehollowsphere(highlightedwithadarkredcolour) issubjected tokinematicboundary conditions,while thevoidsurface(indicated bywhite colour) istraction-free. Inthis

Fig. 4. An illustration of the FE model of the hollow sphere. Symmetry conditions in the three orthogonal directions x 1, x 2, and x 3are imposed. This is indicated by the internal surfaces with a blue colour. The outer surface is subjected to kinematic boundary conditions, while the void surface is traction-free and free to deform. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 5. FE mesh for the hollow sphere model with porosity levels (a) f = 0 . 1 , (b) f = 0 . 01 , and (c) f = 0 . 001 .

work,we haveperformedFEsimulationsforthreedifferentporosity levels(voidvolumefractions) thatcorrespondto f= 0.001,0.01,0.1and the resulting FE models are shown inFig. 5.

Theundeformedconfigurationisusedasthereferenceforthenumericalcomputationstoconformwiththeupper-bound limitanalysisframework,whichisbasedonsmalldeformationtheory.Thus,wedonotdistinguishbetweenthecurrentand referenceconfigurationsinthefollowing.ThevelocityfieldimposedtotheexteriorboundaryoftheRVEreads

u˙ =E˙ ·x (22)

whereE˙ denotesthehomogeneousmacroscopicstrainrateandx=re_r referstoanarbitrarypositionontheouter surface ofthehollowsphere.

Thekinematicboundaryconditionsgoverned byEq.(22)areprescribed suchthat themacroscopicstresstriaxiality(T) andthemacroscopicdeviatoricangle(^{) remainconstant. The}proportional loadingisenforced usingamethodthat was firstproposedbyFaleskogetal.(1998)inaplane-strainsettingandlaterextendedtogeneralthree-dimensionalanalysesby Kimetal.(2004).ArecentpaperbyLiuetal.(2016)elaboratesonthemathematicaldetailsandpresentsageneralunitary transformation matrixthat can beimplemented inan FE model.The formulationadoptedherein is aspecial caseofthe generaltransformationmatrix(Liuetal.,2016)andwaspresentedindetailbyDæhlietal.(2017a).Forthesakeofbrevity, weonlypresentthefinalformofthekinematicconstraintequations.Thesearequitegenerallywrittenontheform

˙ ui=xi

3

k=1

Qik

(

^T,

)

^u˜˙k (23)

where x_i label positionson the outer surface, Q_ik(T,^{) are thecomponents ofan orthogonal} transformation matrixthat explicitlydependsuponthemacroscopicstresstriaxialityandthemacroscopicdeviatoricangle,andu˜_karethedisplacement components ina fictitious system. Note that the Einstein summationconvention isabandoned in Eq.(23). The fictitious systemis inthe sequel referred toas adummy node duetothe useof an externalnode (reference point inAbaqus)to includethesedegrees-of-freedomintheFEmodel.TheconstraintsgiveninEq.(23)relate theglobaldisplacements ofthe unitcellmodeltothedummynode.Theboundaryconditionsarethenprescribedsuchthatu˜˙₁>0,whilethetworemaining displacementsu˜2andu˜3 remainunspecifiedtoaccommodatezeroforcesinthex˜2 andx˜3directions,respectively.

Asasidenote,we mentionthattheMPCusersubroutineisactuallynotneededfortheseparticularanalyses sincethe undeformedgeometryisusedasthereferenceconfiguration.ThekinematicconstraintsgivenbyEq.(23)areconsequently linearandcanbe imposed directlyusinglinear constraintsthroughthe ^∗EQUATIONkeyword inAbaqus.However, a veri- fiedandvalidatedMPC usersubroutinewasalreadyavailableandwasthepreferredoptionforthenumericalsimulations conductedherein.

Themacroscopicstressstate attainedinthelimitanalysisisevaluatedfromtheresultsfiles(.odb) bydirectcalculation ofthevolumeaverage,viz.

= 1 V

Nip

k=1

σ

kV_k (24)

whereVisthetotalvolumeoftheunitcell,V_kand

σ

karethevolumeandthelocalstresstensorofintegrationpointk,and N_ipdenotesthetotalnumberofintegrationpoints.Subsequently,thevonMisesequivalentstressandthehydrostaticstress areevaluatedfrom

^eq=

3

2

^:

h=1

3tr

( )

=

−

h1 (25)

where1isthesecond-orderidentitytensor.

5.2.RemarksontheFEmesh

The three different models corresponding to the porosity levels f=0.001, f=0.01, and f=0.1 were spatially dis- cretizedusing2520,1260,and1050elementsintotal,respectively.Theseconfigurationscorrespondtotheunitcellsshown inFig.5.Theelementsemployedwerequadraticbrickelementswithreducedintegration(C3D20RinAbaqus/Standard).A meshrefinementstudywasinitiallyconductedtoensurethatthelimitloadattainedinthenumericalcalculationshadcon- verged.Tothat end,we examinedthree differentelement types,namelylinearbrickelementswithreduced(C3D8R)and selectively-reduced(C3D8) integration, andquadratic brickelementswithreduced integration(C3D20R).We also limited themesh refinementstudyto an RVEcomprisinga voidvolume fractionof f=0.01 using thequadraticmatrix material (m=2) and astress state corresponding to T=1and=0^◦. Inthe final calculations, we chosea slightlymore refined meshthanwasfoundnecessaryforconvergencefromthismeshrefinementstudy.

Fig.6showstheconvergenceofthe limitloadforvarious spatiallydiscretized unit cellconfigurations.The numberof elementsalong the voidedges and thenumber ofelementsalong the ligament betweenthe voidandthe outer surface (void/ligament)correspondto(i)6/10,(ii)8/12,(iii)10/15,(iv)13/20,and(v)20/30.Thisledtoatotalnumberof260,516, 1020,2560,and9330 elements, respectively.We notethat a seedingbias ratioof10wasused alongtheligament edges, such that the element“height” decreases monotonicallyin the radial directiontowards the voidsurface. This isinferred fromFig.5.WereadilyperceivefromFig.6thatthequadraticelementsgiveconvergedresultsmuchsoonerthanthelinear elements(bothwithreducedandselectively-reducedintegration).Consequently,thequadraticbrickelement(C3D20R)was ourpreferredchoicefortheanalysesconductedherein.

Fig. 6. Convergence rate of the limit load in terms of normalized von Mises equivalent stress ( eq/ σ0) evaluated from a mesh refinement study. The unit cell configuration corresponds to f = 0 . 01 and the loading state to T = 1 and = 0 ^◦in the case of a von Mises matrix ( m = 2 ).

5.3. Imposedmacroscopicstressstates

Wehaveassignedaratherlarge numberofstressstatesinthisworktomap upthenumericalyield surfaces.The im- posed loading conditions correspond to stress triaxiality levels of T=0,1/4,1/2,±1,±2,±10,50 and deviatoric angles of=0^◦,5^◦,10^◦,...,60^◦.Thisyields130differentstress statesforeverycombinationofyieldsurfaceexponentandunit cell configurationandthetotal numberofanalyses conductedforthecurrentwork then amounts1170.Eventhough the analysesare individuallyrathercheapduetothelownumberofelementsandtheexclusionofnonlineargeometryinthe simulations, they collectively representa quite exhaustive computational effort. Also, the use ofa UMAT subroutine and thenon-quadraticmatrixyieldsurfacebothincreasethecomputationaltimefortheanalyses.Wethereforedecidedtoonly include three negativestress triaxiality ratios, whichis deemed sufficient to elucidate the differences betweenthe yield surfacesobtainedfornegativeandpositivehydrostaticstressstates.

6. Resultsanddiscussion 6.1. FEresults

Fig.7showsthenumericalyieldlimitsofthehollowspheremodelwithaporosityof f=0.01intermsofthenormalized macroscopicvonMisesequivalentstress(eq^vm/

σ

0)against themacroscopicdeviatoricangle(^{).Theresultspertaintothe} fullrangeofdeviatoricanglesimposedtotheFEmodel.Resultsforallthethreedifferentyieldsurfaceexponentvaluesare shownandthedatapointsarelabelledbyareddotform=2,agreendiamondform=8,andabluesquareform=20. Thefigurescorrespondto(a)T=−10,(b)T=−1,(c)T=0,(d)T=1,(e)T=2,and(f)T=10.

FromFig.7(c)wereadilyobservethattheshapeoftheyieldsurfacecorrespondstotheunderlyingmatrixyieldsurfacein thecaseofapurelydeviatoricloadingstate.However,theyieldlimitsarereducedbyafactorthatscaleswiththeporosity.

Whenthestresstriaxialityincreases,theyieldstrengthisgreateringeneralizedcompression(=60^◦)thaningeneralized tension(=0^◦). Thisobservationis readilymadefromFig. 7(d)–7(f)andapplies to allthe yield surfaceexponents.The oppositetrendisobservedfornegativestress triaxialities,i.e.theyieldstrengthisgreaterforgeneralizedtensionthan for generalizedcompression,whichisinferredfromFig.7(a)and7(b).Thiseffectcausesasix-foldsymmetryofthemacroscopic yieldsurfacewithrespecttothe stressstate andentailsthatthesymmetry around=30^◦,whichapplies tothematrix yieldsurface,islost.

Fig. 7. Yield points in terms of the normalized von Mises equivalent stress ( ^vmeq/σ0) from the FE analyses of a hollow sphere with porosity f = 0 . 01 . The stress triaxiality ratios correspond to ( a ) T = −10 , ( b ) T = −1 , ( c ) T = 0 , ( d ) T = 1 , ( e ) T = 2 , and ( f ) T = 10 . The results for T = 0 are similar to the underlying matrix yield surface up to a size reduction factor governed by the porosity (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.).

Inthecaseofthequadraticmatrixyieldsurface(m=2),theyieldstrengthoftheRVEattainsaminimumforgeneral- izedtensionforpositivestress triaxialityratiosandincreasesmonotonicallywiththedeviatoricangletowardsgeneralized compression.Theoppositetrendisfoundfornegativestress triaxialityratios.Thus,theexistenceofamaterialheterogene- ityin the formof a void introduces a J₃ dependencein the yield condition forthe porous solid even when the matrix yieldsurfaceonly dependson J₂.Similar results havebeenreported innumericalstudies inthe literature (Cazacu et al., 2013;Keralavarma,2017;Thoré et al.,2011).ThisJ₃ dependence is notpredictedbytheGursonmodel,butcanbecaptured ifhigher-order termsare included in the seriesexpansion of themicroscopic (matrix) dissipation function (Leblondand Morin,2014)orifthemicroscopic dissipationisexactlyintegratedoverthe spatialdomain(Benallal,2017;2018;Benallal etal., 2014; Cazacu etal., 2013). However, we should note that theinfluence of J₃ on the macroscopicyield strength is rathersmallfor thequadratic matrixyield surface (^m=2),andis mostlikely ofsecond-order importance forstructural applications,eventhoughitcanleadtoquitepronounced differencesintermsofporosityevolutionforlargedeformations (Alvesetal.,2014).

Whenthematrixmaterialhasanon-quadraticyieldfunction(m=8orm=20),themacroscopicJ₃ dependenceofthe poroussolid resemblesthat oftheunderlying matrixmaterialforshear-dominated loadingstatesandup tointermediate levelsof stress triaxiality(e.g.|T|∼1).Consequently, theyield strengthis minimizedaround generalizedshear(=30^◦), whichisreadilyseenfromFig.7(b)–7(e).Wenotethatthedifferenceinyieldstrengthbetweengeneralizedcompressionand tensionprevailsfornon-zerostresstriaxialityratios,buttheeffectofJ₃associatedwiththeunderlyingmatrixyieldsurface is more protrusive. However, when the magnitude of the stress triaxiality ratio becomes sufficiently high (e.g. |T|∼10), weobservesimilarbehaviour that wasfoundforthequadraticmatrixyield surface,governedby amonotonic increaseor decreaseinyieldstrengthwiththedeviatoricangledependinguponthesignofthehydrostaticstress(orequivalentlythe stresstriaxiality).Ifwecomparetheyieldstrengthforoppositestresstriaxialityratios,i.e. ±T,theyarerelatedaccordingto eq^vm(^T,)=eq^vm(−T,60^◦−),whichisobservedfromFig.7(a)and7(f)or7(b)and7(d).Thisfeatureofthemacroscopic yieldsurfaceentailscentro-symmetry,forwhichthedefinitionmaybewrittenas

(

eq^vm,

h,

)

=

(

eq^vm,−

h,60^◦−

)

⁽²⁶⁾

Thisis ageneral propertyfor porousplastic materials governed bymatrix yield functionsthat are evenfunctionsof the stress state (Cazacu et al., 2019). From Fig. 7(a) and7(f) it also appears that the difference betweenthe yield strength ingeneralized tensionand generalizedcompression becomesgreater when the exponentm ofthe matrix yieldfunction increases.Itisworthwhiletomentionthateventhoughweobserveayieldstrengthdifferencebetweengeneralizedtension andgeneralizedcompressionforthestresstriaxialityratioswithhighestmagnitude(

|

^T

|

=10),thisdifferenceisonly1%−2%

dependingupon the matrix yield surface exponent.Consequently, there is only a slighteffect of J₃ on the macroscopic yieldingforhighstresstriaxialityratiosevenwhenthematrixmaterialisgovernedbyaJ₃-dependentyieldsurface.

Toshedmorelighton theappearanceofthemacroscopicyield surfacefornon-quadraticmatrixyield surfaces,Figs. 8 and9showplotsoftheyieldlociobtainedfromtheFEunitcellanalysesform=8andm=20,respectively.Thedepicted yieldlocicorrespondtolevelcurves(contours)ofconstantmacroscopicstresstriaxialityprojectedontothedeviatoricplane.

Hence,the datapointson ayield locusaregenerallynotlying onthe samedeviatoricplane intheprincipal stressspace (Keralavarma, 2017). We note that only yield points for=0^◦,5^◦,...60^◦ were evaluated and thesepoints were subse- quentlymirroredtomapuptheentireyieldlocus,utilizingtheisotropyoftheporousplasticsolid.Theresultspresentedin thesefigures correspondtothethree differentporosity levelsusedintheFEanalyses andarelabelledasfollows: f=0.1 (solidbluelines withcirclemarkers), f=0.01 (dotted red lineswithsquare markers), and f=0.001(dashed blacklines withdiamondmarkers).Figs.8(a)–8(i)and9(a)–9(i)pertaintoincreasing valuesofthemacroscopicstresstriaxialityratio:

T=−10,−2,−1,0,0.25,0.5,1,2,10.

Whenthe magnitudeofthe macroscopicstress triaxialityisclose tozero(see Figs. 8(d)–8(f)and9(d)–9(f)),the yield lociareverysimilar inshapetotheunderlyingmatrixyield surface,whichisretainedfora macroscopicstresstriaxiality ratioofT=0.ThisisreadilyseenfromFigs.8(d)and9(d).Thus,inthecaseofpureshearstress states,theyield surface oftheporousductilesolidcoincideswiththeyieldsurfaceofthematrixmaterial,exceptforayielddomainsizereduction dictatedby theporositylevel.Whenthemagnitudeofthestress triaxialityratioincreasesto

|

^T

|

=1and2(seeFigs. 8(b), 8(c),8(g), 8(h),9(b),9(c), 9(g)and9(h)), the highestporosity level (f=0.1)clearlyresultsin amore roundedtriangular shape,whilethetwo lowerporositylevels(f=0.01,0.001)stilldisplay sharpercorners.Thus, forratherhighbutrealistic stresstriaxialitylevels thatarefrequently encounteredinnumericalsimulationsofstructuralcomponents,theeffectofJ₃ onthemacroscopicyieldingofaporousductilesolidisgreatlyinfluencedbytheporositylevel.Byincreasingthemagnitude ofthestresstriaxialityeven further(e.g.T=±10),we readilyseefromFigs. 8(a), 8(i), 9(a)and9(i) thatthemacroscopic yield lociappear almost circular regardlessof the curvature ofthe underlying matrixyield surface (i.e.the value ofm).

Thisobservationpertainstoall threeporositylevelsandimpliesthat underpredominantmacroscopichydrostaticloading, eitherwithpositive ornegative stress triaxiality,the macroscopicyielding is almost unaffected by thedeviatoric loading condition.However,weshouldkeepinmindthatFig.7(a)and7(f)revealedaslightdifferenceintheyieldstrengthbetween generalizedtensionandcompressionforthehigheststresstriaxialitymagnitudes,butthiseffectissosmallthatitishardly visibleontheyieldlocidepictedinFigs.8(a),8(i), 9(a)and9(i).Moreover,inthehydrostaticlimit(i.e.eq^vm=0)theyield locuscoalescestoasinglepointinthedeviatoricplane;thus,itseemsonlynaturalthattheyieldlociforclose-to-hydrostatic stressstates(i.e.^vmeq ≈0and|T|1)becomenearlycircular.

Fig. 8. Discrete yield loci in the principal stress space for m = 8 and constant stress triaxiality ratios projected onto the deviatoric plane. All three porosity levels employed in this work are shown; f = 0 . 1 (blue solid lines - circle markers), f = 0 . 01 (red dotted lines - square markers), and f = 0 . 001 (black dashed lines - diamond markers). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

From the results presentedin Figs. 8and 9,we conclude that the porosity levelplays akey role inhow quicklythe macroscopicyieldsurfaceevolvesfromthehexagonalshapethat isreminiscentofthematrixyieldsurface(seeFig.2),to the roundedtriangular shape,andfinally tothe nearlycircularshape prevalentforpredominant hydrostaticstress states.

Asimilar finding wasreportedin theprevious studybyRevil-BaudardandCazacu (2014a) foraTrescamatrix behaviour, whichsharesmanyofthesamefeaturesasthematrixyieldsurfacewithm=20usedinthisstudy.Wenotethatthereare rathersmalldifferencesbetweenthemacroscopicyieldpointsfor f=0.01and f=0.001uptoratherhighstresstriaxiality ratios.Roughlyspeaking,we observealmost nodifference betweenyieldlocifor f=0.001and f=0.01 withinarangeof stress triaxialityratioscorresponding to−1≤T≤1(see Figs.8(c)–8(g)and9(c)–9(g)).Thisisactually quitean important rangeofstresstriaxialityratiosandporositylevels intermsofpractical applicationsforstructuralmetal alloys.Theeffect of theporosity becomes greater whenthe magnitudeof the stress triaxialityratio increases.From the numericalresults shownin Figs.8 and9,we observepronounced differencesbetweentheyield lociforall threeporosity levels whenthe magnitudeof the stress triaxialityratio islarger than roughly |T|≈2.In the caseofthe highestmagnitude ofthe stress triaxialityratio(

|

^T

|

=10),theyieldlimitsareverydifferentdependingupontheporositylevel.Assuch,thehydrostaticlimit

Fig. 9. Discrete yield loci in the principal stress space for m = 20 and constant stress triaxiality ratios projected onto the deviatoric plane. All three porosity levels employed in this work are shown; f = 0 . 1 (blue solid lines - circle markers), f = 0 . 01 (red dotted lines - square markers), and f = 0 . 001 (black dashed lines - diamond markers). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

(i.e.eq^vm=0) isgreatlyaffectedby theporositylevel,whichisinaccordancetotheporousplasticitymodelspresentedin Section4.

Wealsofindthattheyieldlocusshapetransitionisinfluencedbythemagnitudeoftheyieldsurfaceexponent(m).This isto some extent observed by comparing the yield lociin e.g. Figs. 8(g) and9(g). However, the influence ofthe matrix yieldsurfaceexponentontheshapetransitionseemstobeonlysecondordertothatoftheporositylevel,andtheporosity isclearlythedominantfactorforthe yieldsurfaceevolution withincreasingmagnitudeof thestress triaxiality.Thiswas alsoaddressedinthestudybyRevil-BaudardandCazacu(2014a)showing3Dyieldsurfacesforconstantlevelsofthemean strain.

Beforewe proceed totheassessment ofthetwo porousplasticitymodels,we wouldlike tohighlighta particularob- servationfromthenumericalanalyses withthequadraticmatrixyieldfunction(m=2).WhenT=0,theporousplasticity modelsdisplayayieldsurfacethatisidenticaltotheunderlyingmatrix,exceptthattheyieldstrengthisreducedbyafactor thatscaleswiththeporositylevel.WithreferencetoFig.7(c),thissuggeststhatweshouldobtainastraightlineform=2. However,fromthenumericalresultsshowninFig.10,wemaynoticethatwedonotobtainexactlyastraightline,although

Fig. 10. Yield points for m = 2 and T = 0 plotted against the macroscopic deviatoric angle. The data points correspond to the von Mises equivalent stress at yielding ( eq^vm) normalized by the von Mises equivalent stress at yielding in generalized tension ( ^vmeq,0= eq^vm(= 0 ^◦)) for all three levels of porosity.

thedeviationisonlyslight.Theyieldstrengthratherdecreasesforgeneralizedshearandmonotonicallyincreasestowards theaxisymmetriclimits.Theobservedeffectisconsistentandpertainstoallthreelevelsofporosity,butitdecreaseswith decreasingporosity.Assuch,wedonotbelievethistobecausedbynumericalerrors,butrathertoresultfromthespatial heterogeneity ofthemechanicalfieldsandthatthenumericalanalyses accommodateamoregeneralsetofvelocityfields comparedtotheuniquetrialvelocityfieldusedintheupper-boundlimitanalysis.

6.2. Assessmentoftheporousplasticitymodels

Inthissection,wecomparethetwoanalyticalmodelspresentedinSections4.1and4.2tothenumericaldataobtained fromtheFElimitanalysesanddiscusstheir respectiveperformance inlight oftheyieldlimitsprovided bytheFEcalcula- tions.The latterareconsideredto betheexactyield limitsforthehollowspheremodel,whilethetwo analyticalmodels areupper-boundsolutions.However,onlythemodeldevelopedbyBenallal(2017,2018)can bea prioriconsidered as a rig- orousupper-boundsolution(seeSection4.1),sincetheheuristicextensionisonlymathematicallyrigorousforthequadratic matrixyieldsurface;i.e.whenitreducestotheoriginalGursonmodel(Gurson,1977).

Figs.11–13compareyieldlocifortheporousplasticitymodelstotheFElimitanalysesforallthreematrixyieldsurface shapesandforallporosity levels examinedherein. Thesolid anddashed linescorrespond totherigorousmodelandthe heuristic model, respectively, while the circle markers indicate the FE yield points. The porosity levels are indicated by different colours: f=0.001 (green), f=0.01 (red), and f=0.1 (blue).We must emphasize that the yield locirepresent contoursofconstantstresstriaxialityandnotlevelsofconstanthydrostaticstress.Assuch,theyieldlociarecurvedsections of theyield surfacein theprincipal stress spaceand theloci derived fromthe FEanalyses, the rigorousmodel,andthe heuristic model generally do not lie in the same plane. However, they correspond to sectionsof the same macroscopic stress triaxiality, which is an importantparameter that is frequently used asa referenceparameter in FEsimulations of materialtestsandstructuralcomponentsandisofkeyimportancefortheporosityevolutioninlarge-straintheory.

From Figs.11–13,we observethattheresultsobtainedwiththerigorousmodelareingoodagreementtothose ofthe FEsimulationsforallstressstatesandshapesofthematrixyieldsurface.Thepredictionsarelessaccurateforhighporosity levels, which conforms with the results of the previous study by Cazacu etal. (2014b), but the shape of the yield loci are stillinvery goodagreementtotheFEsimulations. Weobservethat therigorousmodelreflectsthetransformationof theyieldsurfacewithstresstriaxiality,fromthehexagonalshape (

π

/6-symmetry)atlowstress triaxialityto therounded triangularshape(

π

/3-symmetry)atfairlyhighstresstriaxiality.Thisisperhapsmosteasilyverified fromFig.13.Boththe rigorousmodel andtheheuristic modelindeedprovide upper-boundsolutions forthematrix yieldsurfaces andporosity levelsstudiedherein.Thisisapriorisatisfiedfortherigorousmodelduetotheexactsolutionofthekinematiclimitanalysis, butcannotbeguaranteedfortheheuristicmodelotherthanform=2.

Whiletherigorousmodelgenerallyprovidesratheraccurate predictions,weobservethatalsotheheuristicmodelpre- dicts theyieldstress withgoodaccuracyforall stress statesandshapesofthe matrixyieldsurfacewhenthe porosityis low(f=0.001).Evenfortheintermediateporositylevel(f=0.01)thepredictionsoftheheuristicmodelarerathersimilar tothoseoftherigorousmodelandtheFEresults,althoughsomeminordiscrepanciesarevisible,especiallyforstressstates closetogeneralizedtensionorgeneralizedcompression(seee.g.Figs.12(e)and13(e)).Whentheporosityisveryhigh(i.e.

f=0.1),theaccuracyoftheheuristicmodeldeterioratesforthenon-quadraticmatrixyieldsurfaces.Thelargestdiscrepancy isobservedaroundaxisymmetricstressstates(=0^◦ or60^◦)andintermediatetohighstresstriaxialityratios(

|

^T

|

≈1−2

Fig. 11. Comparison between the unit cell results (solid points), the rigorous model (solid lines), and the heuristic model (dashed lines) for yield surface exponent m = 2 . The plots correspond to porosity levels f = 0 . 001 (green curves), f = 0 . 01 (red curves), and f = 0 . 1 (blue curves) and stress triaxiality levels of (a) T = −10 , (b) T = −1 , (c) T = 0 , (d) T = 1 , (e) T = 2 , (f) T = 10 . (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

intheseanalyses),wheretheheuristicmodelmarkedlyoverestimatestheactualyieldstrengthevaluatedfromtheFEcalcu- lations.Themainreasonforthesedeviationsisthattheheuristicmodelpredictsamacroscopicyieldsurfacethatretainsthe shapeoftheunderlyingmatrixyieldsurfaceregardlessofthestresstriaxialityratioandtheporositylevel.Consequently,the heuristicmodelcannotaccountforthetransformationofthemacroscopicyieldsurfacethatdrasticallychangestheshapeof theyieldloci;i.e.fromthehexagonalshapetotheroundedtriangularshapeforincreasingtriaxiality.Thisisneitherinac- cordancetothepredictionsobtainedintheFElimitanalysesnortothepredictionsoftherigorousporousplasticitymodel, whichbothdisplaythistypeofyieldsurfacetransformation.Sincethistransformationwasfoundtobecloselyrelatedtothe porositylevel(seee.g.Fig.9(g)andthediscussioninSection6.1),whereahighporositylevelfacilitatesearlierdevelopment oftheroundedtriangularshape,thepredictionsoftheheuristic modelareconsequentlyhamperedby increasingporosity level.However,whenthemagnitudeofthestresstriaxialitybecomessufficientlyhigh(e.g.|T|∼10inourresults),thepre- dictionsoftheheuristicmodelisseentoagreewellwithboththeFEyieldpointsandthepredictionsoftherigorousmodel.

Whetherthisisofpracticalrelevanceisofcoursedebatable,sincesuchhighstresstriaxialityratiosarerarelyencountered instructuralapplicationsduetobluntingofthehighlyconstrainedregionsduringdeformation.

Figs. 11–13provide a reasonablecomparison ofthe porous plasticitymodels to the FE dataand clearlyhighlight the predictive capabilities ofthe two models. However, the quality of the predictions is perhaps more easily assessed from Figs.14and15,whichshowhowtheyieldstressintermsofnormalizedvon Misesequivalentstress(eq^vm/

σ

0)varieswith themacroscopicdeviatoricangle(^).

Fig.14(a)–14(f) pertain to all three matrix yield surface shapesat different levels of stress triaxiality. The solid lines correspondtotherigorousmodel,thedashedlinestotheheuristic modelandthepointsindicatedby markerspertainto theFEyieldpoints.Thedifferentmatrixyieldsurfaceshapesarelabelledasfollows:(i)m=2ishighlightedwithredcolour andsolidmarkers,(ii)m=8ishighlighted withgreencolouranddiamondmarkers,while(iii)m=20ishighlightedwith