framework of gradient crystal plasticity Modeling of surface effects in crystalline materials within the Journal of the Mechanics and Physics of Solids

ContentslistsavailableatScienceDirect

Journal of the Mechanics and Physics of Solids

journalhomepage:www.elsevier.com/locate/jmps

Modeling of surface effects in crystalline materials within the framework of gradient crystal plasticity

Xiang-Long Peng

^a,b

, Edgar Husser

^c

, Gan-Yun Huang

^a,∗

, Swantje Bargmann

^b,∗

aDepartment of Mechanics, School of Mechanical Engineering, Tianjin University, Tianjin 300350, PR China

bChair of Solid Mechanics, School of Mechanical Engineering and Safety Engineering, University of Wuppertal, Wuppertal 42119, Germany

cInstitute of Continuum Mechanics and Material Mechanics, Hamburg University of Technology, Hamburg 21073, Germany

a rt i c l e i n f o

Article history:

Received 20 September 2017 Revised 5 December 2017 Accepted 12 January 2018 Available online 12 January 2018 Keywords:

Strain gradient Crystal plasticity Size effect Surface effects Surface yielding

a b s t ra c t

Afinite-deformationgradient crystalplasticitytheoryisdeveloped,whichtakesintoac- count theinteractionbetweendislocations andsurfaces. Themodelcaptures bothener- geticand dissipativeeffectsfor surfacespenetrableby dislocations.By takingadvantage oftheprincipleofvirtualpower,thesurfacemicroscopicboundaryequationsareobtained naturally.Surfaceequationsgovernsurfaceyieldingandhardening.Athinfilmundershear deformation servesas abenchmarkproblemforvalidation oftheproposed model.It is foundthatbothenergeticanddissipativesurfaceeffectssignificantlyaffecttheplasticbe- havior.

© 2018TheAuthors.PublishedbyElsevierLtd.

ThisisanopenaccessarticleundertheCCBYlicense.

(http://creativecommons.org/licenses/by/4.0/)

1. Introduction

In the last two to three decades, many efforts have been devoted to studying the mechanical properties of small- scaledand/orfine-grainedcrystallinematerialswidelyusedinsmall-scaledengineeringsuchasmicroelectronicsandmicro- electromechanicalsystemstoensuretheirperformanceandreliabilityinpracticalapplications.Ithasbeenfoundinvarious experiments(seee.g.,Flecketal.,1994;Guetal.,2012;McElhaneyetal.,1998;Swadeneretal.,2002)thatcrystallinema- terialsatsmallscalesusuallypossessplasticbehaviorsquitedifferentfromthoseoftheirbulkcounterparts,amongwhichis thesize-dependenceoftheyieldstressandstrain hardeningbehavior.Sinceconventionalplasticitytheoriesareessentially size-independent,thesesizeeffectsarebeyondtheir predictability.Asiswell known,themacroscopicplasticdeformation incrystalsmainlyresultsfromthemicroscopicdislocationglidingontheindividualcrystallographicslipplanes,thusthose unusualplastic behaviorsare closelyrelevantto thedistinctive dislocationactivities insmall-scaled crystals.Accordingto Ashby (1970),dislocations incrystals can be classifiedinto two types,namely statisticallystoreddislocations (SSDs) and geometricallynecessarydislocations(GNDs).Thelatterisstoredincrystalstoaccommodatelatticecurvature.Asdiscussed byNixandGao(1998),GNDsmayberesponsibleforsizeeffectsduetostraingradientssincetheyaremorelikelytoaccu- mulatewheninhomogeneousplasticdeformationtakesplace.Furthermore,withthedecreaseofgrainsizeand/orgeometry sizeofcrystals,thevolumetograin boundary(GB)/surfaceratioincreasesdramatically,whichmaypromptdislocationsto reachandtheninteractwiththeGB/surface.Forinstance,asdemonstratedbyatomisticsimulations(seee.g.,Spearotetal.,

∗ Corresponding authors.

E-mail addresses: gyhuang@tju.edu.cn (G.-Y. Huang), bargmann@uni-wuppertal.de (S. Bargmann).

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

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

( http://creativecommons.org/licenses/by/4.0/ )

2005;VanSwygenhovenetal.,2002;Yuasaetal.,2010),grainboundariesactassourcesandsinks,emittingandabsorbing dislocations.Insinglecrystals,dislocationsareobservedtoescapetowardsandnucleateatthesurface,whichdramatically changesthestrainhardeningbehavior(Ohetal.,2009;Shanetal.,2008;Wangetal.,2012).Inordertoconstructpredictive theoriesfortheplasticbehaviorinsmall-scaledcrystallinematerials,thedislocation-relevantmechanismsbothinthebulk andattheGB/surfaceneedtobeproperlyconsideredonthecontinuumlevel.

Sofar,sincethepioneeringworkbyAifantis(1984),numerousstraingradientplasticitytheorieshavebeenproposedby incorporatingtheplasticstraingradientortheGNDdensity,andhavesuccessfullypredictedexperimentallyobservedsize- effects(seee.g.,ChenandWang,2000;Flecketal.,1994).AsdiscussedinKurodaandTvergaard(2008),theexistingstrain gradienttheoriescanbeclassifiedintotwotypesnamelythelower-ordertheories(seee.g.,AcharyaandBassani,2000;Chen andWang,2000;Eversetal., 2002)andthehigher-ordertheories(seee.g.,Bargmannetal.,2011;Evers, 2004;Fleckand Hutchinson,1993;1997;Gromaetal.,2003;Gurtin,2000;2002;Yefimovetal.,2004)basedonhowtheeffectsofGNDsare considered.Ofparticularinterestarehigher-ordertheories,inwhichinadditiontothetraditionalforcebalanceequationand thetractionboundarycondition,themicroscopicgoverningequationsandthecorrespondingmicroscopicboundarycondi- tionsareinvolved.¹Thehigher-ordertheoriesprovidethepossibilitytomodeldislocationactivitiesataGB/surfacethrough properlychoosingmicroscopicboundaryconditions.Forsimplicity,themicrohardandmicrofreeconditionsrepresentingtwo extremecasesareusuallyimposed,andtheunderlyingassumptionisthattheGB/surfaceiscompletelyimpenetrableresp.

freelypenetrabletodislocations.However,thesetwoextremeGB/surfacemodelsarenotsufficienttocapturedislocationac- tivitiesattheGB/surface.Todate,afewintermediateGB/surfacemodelshavebeenproposed.CermelliandGurtin(2002)de- velopedadissipativegrainboundarymodelwhereavisco-plasticrelationforgrainboundarymicrotractionisassumed.By consideringtheinterfaceenergycontributedbytheinteractionbetweendislocationsandgrainboundaries,otherresearchers proposedenergeticgrainboundarymodels(seee.g.,Aifantisetal.,2006;FredrikssonandGudmundson,2005;Gudmundson, 2004).Thesemodelsarepurelyphenomenological.Gurtin(2008)constructedagrainboundarymodelwheretheeffectsof grainmisorientationandgrainboundaryorientationareconsideredthroughtheexactlyderivedgrainboundaryBurgersten- sor.InvanBeersetal.(2013),agrainboundarymodelsimilartothatofGurtin(2008)wasputforward,inwhichavectorial measureofthedensityofgrainboundarydefectwasintroduced.AsillustratedbyGottschalketal.(2016),thesetwomod- elsareessentiallyidenticalifconsideringplanarproblems.vanBeersetal.(2015)extendedthepreviousmodel(vanBeers etal.,2013) byfurtheraccountingfortheredistributionofthegrainboundarydefects.Ekhetal.(2011)suggestedamicro- flexibleboundaryconditionbywhichthecrystallographicmisorientationbetweenadjacentgrainsistakenintoaccount.The thermaleffectsonthemodelarefurtherconsideredby BargmannandEkh(2013).Inthegrainboundarymodelproposed byWulfinghoff etal.(2013),acriterionforinterfaceyieldingisintroducedtoaccountfortheresistanceofgrainboundaries againstslipoccurrence.Kuroda(2017)proposedastrategyformodelingvariousmicroscopicinterfaceboundaryconditions byintroducingascalarquantitytocontroltheslippingrateattheinterface,basedonwhichanorientation-dependentgrain boundarymodelisproposedandthenimplementednumericallytostudytheplasticbehaviorofbicrystalmicropillars.

Somemodelsstudydislocation-surfaceinteractionsinsinglecrystals.HuangandSvendsen(2010)proposedanenergetic surfacemodelbyconsideringthechangeofsurfaceenergycontributedbysurfacestepsduetodislocationabsorptionbysur- faces.Similarly,HurtadoandOrtiz(2012)derivedanexplicitexpressionforthedensityofsurfaceenergybyconsideringthe relativeorientationbetweenslipsystemsandsurfaces.Inbothmodels,thedensityofsurfaceenergydependslinearlyonthe plasticslipsatthesurface.Byfurtherconsideringtheinteractionenergybetweensurfacesteps,PengandHuang(2015)con- structedaphysicallybasedenergeticsurfacemodelwithintheframeworkofwork-conjugatedhigher-ordergradientcrystal plasticity.Inthismodel,asurfaceyieldingconditionisderivednaturally,andinteractionsbetweenslipsystemsareconsid- ered.Recently, Husseretal.(2017) suggestedaflexible surfaceboundarymodelby establishinga directrelationbetween thesurfacesliprateandthemicroforcevectorprojectedontotheboundary.Usingasurfaceyieldcriterion,anon-idealized surfacebehaviorwasaddressedincludingdissipativeeffectsofsize-dependentsurfaceyielding.

Asreviewedabove,modelingGB/surfaceeffectshasrecentlybeentheforefrontinthefieldofcrystalplasticityatsmall scales.Yet,intheexisting GB/surfacemodels,lessattentionhasbeenpaidontheunderlyingphysicalmechanisms.Toin- corporatetheunderlyingphysicalmechanisms relevanttointeractions betweendislocationsandsurfacesasillustratedby insituexperiments intothecontinuum model, werevisit modeling dislocation absorption bysurfaces.On theone hand, afterdislocationabsorptionbysurfaces,surfacestepsform(seee.g.,Shanetal.,2008;Wangetal.,2012;Zhongetal.,2017), whichresultsinthechangeofsurfaceenergy. Thus,suchsurface-drivenchangeinthefree energymustbeproperlycon- sidered.Ontheotherhand,asisobservedinOhetal.(2009),dislocationsstopforawhilenearasurfacecontaminatedby anoxidelayerbefore theyareabsorbed.Thisisstrongevidenceforthe existenceofanincreasedresistanceagainstdislo- cationabsorptionatsurfacesaffectedbyoxidelayers,coatingsorinitialdefects.Furthermore,accordingtotheexperiment byZhongetal.(2017),dislocationabsorptionbysurfacesmaybeviewedasadissipativeprocess.Allthoseindicatethedis- sipativenatureofdislocationabsorptionby surfaces.Inaddition,modelsconfinedtosmalldeformationareinadequatefor situationswithsevereplasticdeformation.Consequently,anadvancedsurfacemodelwithenergeticanddissipativesurface effectsisproposed intheframework offinitedeformationgradientcrystalplasticity.Energeticsurfaceeffectsaremodeled byextendingthephysicallymotivatedsurfacemodelproposedinPengandHuang(2015).Toconsiderdissipativesurfaceef- fects,arate-dependentconstitutiverelationforthedissipativesurfacemicroforceisproposedwithinthenewsurfacemodel.

1For a unified formulation and comparison of higher-order theories, see ( Svendsen and Bargmann, 2010 ).

Thetheoryisnumericallyimplementedbyusinganin-housefiniteelementcodebasedonadual-mixedfiniteelementso- lutionstrategytostudyabenchmarkproblem.

2. Modelframework

2.1. Basickinematicsoffinitedeformationcrystalplasticity

Infinitedeformationcrystalplasticity,thebasicassumptionisthatthedeformationgradienttensorFismultiplicatively decomposedintoanelasticpartF^eandaplasticpartF^p asfollows

F=F^e·F^p. (1)

Theplastic partF^p describestheplasticdeformation fromthereferenceconfigurationtothestress-freeintermediatecon- figuration resultingfromdislocation glidingon thecrystallographic planes.Theelastic partF^e by whichthe intermediate configurationisdeformedtothecurrentconfigurationcomprisesthelatticedistortionandrotation.Itisassumedthatthe plasticdeformationdoesnotchangethevolumesuchthatdetF^p=1withdetdenotingthedeterminantofatensor.Inthe following,ifnecessary,thesubscriptsiandcareusedtoidentifyquantitiesintheintermediateandthecurrentconfigura- tion,respectively.

Basedonrelation(1),thevelocitygradienttensorLisdecomposedas

L=

∇

^{cu˙ =F˙·F−1=Le+Fe·Lp·Fe−1, (2)}

withtheelasticvelocitygradientL^eandtheplasticvelocitygradientL^p definedas

L^e=F˙^e·F^e−1 (3)

and

L^p=F˙^p·F^p−1, (4)

where[

∇

^a]i j=a_i,jdenotesthegradientofavectora,thesuperposeddotdenotesthematerial-timederivative,thesuper- script −1 denotesthe inverseof a tensor,andu isthe displacement vector. Letusconsider a crystalinwhich each slip system

α

^{isidentifiedbyaslipdirections(α)andaslipplanenormalm(α)attachedtothelatticespaceinthe}intermediate configuration.Theisoclinicassumption isadoptedsuchthat s^(α) andm^(α) remainunalteredfromthereferenceconfigura- tiontotheintermediateconfiguration.Then,bydenotingthesliprateintheintermediateconfigurationas

γ

˙^{(α),theplastic} velocitygradienttensorisexpressedasthesumoverthecontributionsfromallslipsystems

L^p= α

γ

˙^(α)s(α)m^(α). (5)

Therelevantstrain measures,namelytheGreen–Lagrangestrain tensor,andtherightCauchy–Green stretchtensorand theirelasticpartsarerespectivelydefinedas

E=1

2[C−I], C=F^T·F (6)

and

E^e=1

2[C^e−I], C^e=F^eT·F^e, (7)

wherethesuperscriptTdenotesthetransposeofatensorandIistheidentitytensor.

TheevolutionofedgeandscrewGNDdensities

ρ

_i^ge(α)and

ρ

_i^{gs(α)areputdownas}

ρ

˙_i^ge(α)=−1

b

∇

i

γ

˙^(α)·s^(α) (8)

and

ρ

˙i^gs(α)=1

b

∇

i

γ

˙^(α)·p^(α), (9)

whereb is themagnitudeofthe Burgersvector, andp^(α)=s^(α)×m^(α).Theinitial conditionsare

ρ

i^ge(α)

|

t=0=

ρ

0i^ge(α) and

ρ

i^gs(α)

|

^t=0=

ρ

0i^gs(α)with

ρ

^ge0i^(α)and

ρ

0i^gs(α)beingtheinitialGNDdensities.

2.2. Balanceequationsandboundaryconditions

Inthefollowing,themacroscopicandthemicroscopicbalanceequationsarederived viatheprincipleofvirtualpower.

FollowingGurtin(2002),intheintermediateconfiguration,theinternalpowerinthebulkisassumedtobeexpendedbythe second Piola–Kirchhofstress tensorS^e power-conjugatetotherateoftheGreen–Lagrange straintensorE˙^e,themicroforce

π

i^(α) power-conjugate to the slip rate ˙

γ

^(α) and the microstress

ξ

i^(α) power-conjugateto theslip rategradient

∇

i

γ

˙^{(α) for}

eachslipsystem

α

^{.Inthepresentwork,anadditional}contributiontotheinternalpowerexpendedbythemicroforce

η

^(α)_i

power-conjugate to the sliprate

γ

˙^{(α) at the surfaceS}i is introduced to consider the effect of dislocation absorption by surfaces.Thus,thetotalinternalpowerP^intisexpressedintheintermediateconfigurationas

P^int=

Vi

S^e:E˙^edVi+ α

Vi

π

i^(α)

γ

˙^(α)+

ξ

^(α)i ·

∇

ⁱ

γ

^˙(α)

dVi+ α

Si

η

i^(α)

γ

˙^{(α)dSi, (10)}

wherethesecond Piola–Kirchholfstress tensorS^e isdefinedasS^e=F^e−1·J

σ

·F^e−T with

σ

^{beingthe Cauchystress tensor.}

Accordingly,theexternalpowerP^extintheabsenceofbodyforcesiswrittenas P^ext=

Si

T_i·u˙dS_i+ α

Si

^(α)i

γ

˙^(α)dSi, (11)

whereT_iistheprescribedtractionand^(α)_i ^isthemicrotractiondefinedintheintermediateconfiguration.Accordingtothe principleofvirtualpower,the variationoftheinternalpowerwithrespecttothevelocityu˙ andthesliprate

γ

˙^{(α) equals} thatoftheexternalpower,whichgives

Vi

S^e:

δ

^E˙edVi+ α

Vi

π

i^(α)

δ γ

^˙(α)+

ξ

^(α)i ·

∇

i

δ γ

^˙(α)

dV_i

+

α

Si

η

i^(α)

δ γ

^˙(α)dSi=

Si

T_i·

δ

^u˙dSi+ α

Si

^(α)i

δ γ

^˙(α)dSi. (12) FromEqs.(2),(3),(5)and(7),thevariationofE˙^eisexpressedas

δ

^E˙e=1₂

F^eT·

∇

i

δ

^{u˙ −}

α

C^e·s^(α)m^(α)

δ γ

^˙(α)

+1 2

F^eT·

∇

i

δ

^{u˙ −}

α

C^e·s^(α)m^(α)

δ γ

^˙(α)

T

. (13)

Then,by substituting Eq.(13)intoEq.(12) andusingthe divergencetheorem,theprinciple ofvirtual power isrewritten

as

Vi

Div_i(F^e·S^e)·

δ

^u˙dVi− α

Vi

Div_i

ξ

⁽i^α)+s^(α)·M·m^(α)−

π

i^(α)

δ γ

^˙(α)dVi

+

α

Si

ξ

i^(α)·N_i+

η

⁽i^α)−

⁽i^α)

δ γ

^˙(α)dSi+

Si

[[F^e·S^e]·N_i−T_i]·

δ

^u˙dSi=0, (14)

whereN_i is the surfaceoutward normal vector, Div denotes thedivergence operator, andthe Mandelstress tensor M is expressedas

M=C^e·S^e. (15)

First, considering the special case without plastic slip (

δ γ

^˙(α)=0), the requirementthat Eq. (14)should be satisfied for arbitrary

δ

^{u˙ results in the} macroscopicforce balance equation andthe corresponding standard traction conditionin the intermediateconfiguration,i.e.,

0=Div_i

(

^Fe·Se

)

⁽¹⁶⁾

and

[F^e·S^e]·Ni=Ti (17)

atthepartofsurfaceS^T_i withprescribedtraction.Then,giventhevalidityofEq.(14)forarbitrary

δ γ

^{˙(α)inthecasewithout}

macroscopicdisplacement (

δ

^u˙ =0), the microscopicbalance equation andtheassociated microscopic boundarycondition foreachslipsystem

α

^{areobtainedasfollows}

Div_i

ξ

i^(α)+

τ

i^(α)−

π

i^(α)=0, (18)

with

τ

i^(α)=s^(α)·M·m^(α)beingtheSchmidstress,and

ξ

i^(α)·N_i+

η

i^(α)−

^(α)_i =0. (19)

2.3.Constitutiverelations:theorywithsurfaceeffects

Theelasticbehavior ofthe crystalisassumedtobe compressibleneo-Hookeansuchthat thestrainenergydensity

ψ

i^e

hasthefollowingform

ψ

i^e

(

^Ee

)

⁼

μ

^Ee:I+

λ

2[ln

(

^J

)

^]2−

μ

^ln

(

^J

)

⁽²⁰⁾

withJ=

det(^2Ee+I)^,and

λ

^and

μ

^beingtheLamé constants.Then,thehyperelasticconstitutiverelationforS^eisderived as

S^e=

∂ψ

i^e

(

^Ee

)

∂

^{Ee =}

μ

^I+[

λ

^ln

(

^J

)

⁻

μ

^{]Ce−1. (21)}

Generally,the microstress

ξ

^(α)_i ^{maycomprise anenergetic partanda} dissipativepart(Bargmannet al.,2014;Gurtin and Anand, 2005).Here, forsimplicity,

ξ

_i^{(α) isassumedtobepurely energetic.AccordingtoKuroda andTvergaard(2008)and}

Ertürk etal. (2009),the energetic microstress

ξ

^(α)i inthe work-conjugate theories can be essentially relatedto the back stress

τ

_i^b(α)inthenon-work-conjugatetheories(seee.g.,Bayleyetal.,2006)asfollows

τ

_i^b(α)=−Divi

ξ

_i^(α). (22)

In Ertürk et al. (2009), based on relation (22) and the expression of the back stress

τ

_i^{b(α) derived by}

Bayley et al. (2006) through considering the elastic interactions between GNDs within the same slipsystem and those fromdifferentslipsystems,theconstitutiverelationfor

ξ

i^(α)isgiven.Here,wedirectlyadopttherelationtherein,thus

ξ

i^(α)

isexpressedas

ξ

i^(α)=

μ

^bl2

8[1−

ν

^]

β

ρ

i^ge(β)B^e(αβ)+

μ

^bl2

4

β

ρ

i^gs(β)B^s(αβ) (23)

with

B^e(αβ) =

3

s^(α)·s^(β)

m^(α)·s^(β)

+

s^(α)·m^(β)

m^(α)·m^(β)

+4

ν

s^(α)·p^(β)

m^(α)·p^(β)

m^(β)

−

s^(α)·s^(β)

m^(α)·m^(β)

+

s^(α)·m^(β)

m^(α)·s^(β)

s^(β), B^s(αβ) =−

s^(α)·p^(β)

m^(α)·s^(β)

+

s^(α)·s^(β)

m^(α)·p^(β)

m^(β) +

s^(α)·s^(β)

m^(α)·m^(β)

+

s^(α)·m^(β)

m^(α)·s^(β)

p^(β), (24)

where

ν

^{isPoisson’s ratioandldenotes theradiusof thedomainwithin whichthe}interaction betweenGNDsis consid- ered.Alternatively,theinteractionbetweendifferentslipsystemscanbeaccountedforintherelationforGNDdensitiesin Eqs.(8)and(9)asdoneinBargmannetal.(2011)andHusser&Bargmann(2017),butthiswayisnotpursuedhere.

Thescalarmicroforce

π

_i^(α)beingdissipativeinnatureisassumedtofollowavisco-plasticpower-law,i.e.,

π

i^(α)=R^b_i^(α)

γ

˙^(α)

γ

˙0^b

mb

sgn

γ

˙^(α)

, (25)

wheresgn()denotesthesignfunction,m_b istherate-sensitivityexponentinthebulk,and

γ

˙₀^{bisthereferencesliprate.The} slipresistanceR^b_i^(α)resultingfromtherandomtrappingprocessesbetweenSSDshasthefollowinglinearhardeningform,

R^b_i^(α)=R^b₀^(α)+ β

H^b(αβ)

γ

acc^(β), (26)

where

γ

acc^(β)=

γ

˙^(β)

dt is the accumulated plastic slip. The initial slip resistance R^b₀^(α) may depend on the initial SSD density,andH^b(αβ) isthelocalhardeningmodulus.BysubstitutingEq.(25)intoEq.(18),themicroscopicbalanceequation isrewrittenastheevolutionequationofplasticsliprate,i.e.,

γ

˙^(α)=

γ

^˙0^b

Divi

ξ

_i^(α)+

τ

_i^(α)

R^b_i^(α)

_m¹

b

sgn Div_i

ξ

i^(α)+

τ

i^(α)

, (27)

The microforce

η

_i^{(α) atthe surfaceisdivided intoan energetic part}

η

^en_i ^{(α) anda} dissipativepart

η

^dis_i ^{(α).According to}

Pengand Huang(2015),the powerexpended by the energeticmicroforce

η

^en_i ^{(α) atthesurface equals thechange ofthe}

surfaceenergyresultingfromsurfacestepsformedafterdislocationabsorption,

α

Si

η

i^en(α)

γ

˙^(α)dSi−

Si

ψ

˙i^SdS_i=0, (28)

where

ψ

˙i^S is the rate of the surface energy densityresulting from surface steps formed after dislocation absorption. In PengandHuang(2015),theexpressionforthedensityofsurfaceenergyinthesmalldeformationcasewasstrictlyderived by considering both the surfaceenergy change due to the increase of surfacearea andthe interaction energy between surfacesteps. Following the ideatherein andconsidering the effectof finitedeformations, theexpression for

ψ

˙i^S can be readilyderived,andwithoutgoingintothedetails,theresultisdirectlygivenhere,i.e.,

ψ

˙i^S=

α

1^(α)

γ

˙^(α)

s^(α)·N_i

sgn

γ

^(α)

+ α,β

2^(αβ)

s^(α)·N_i

m^(α)×N_i

s^(β)·N_i

m^(β)×N_i

γ

˙^(α)

γ

^{(β). (29)}

Thefirstandsecondtermsrepresentcontributionsfromtheincreaseofsurfaceareaandtheinteractionbetweensurface steps,respectively.^(α)1 denotes the surface freeenergy per unit surface step area and^(αβ)2 is the surface hardening mod- ulusmeasuringtheinteractionstrengthbetweensurfacestepsfromslipsystems

α

^and

β

^{.InEq.(29),thesurfaceoutward}

normalvectorN_iintheintermediateconfigurationdependsondeformation,and,hence,theexplicitexpressionfortheden- sityofsurfaceenergyisnotavailable,whichisdifferentfromthesmalldeformationcasewherethesurfaceoutwardnormal vectorremainsunaltered.Then,fromEq.(28),

η

i^en(α)isexpressedas

η

_i^en(α)=

^(α)1

s^(α)·Ni

sgn

γ

^(α)

+ β

2^(αβ)

s^(α)·Ni

m^(α)×Ni

s^(β)·Ni

m^(β)×Ni

γ

^(β). (30)

Thedissipativesurfacemicroforce

η

^disi ^(α)isintroducedtoconsiderthedissipativenatureofdislocationabsorptionbysur- facesasindicatedbyinsituexperiments(Ohetal.,2009;Zhongetal.,2017).Toconstructaconstitutiverelationfor

η

^dis_i ^(α)

physically,detailedinformationontheprocess fora dislocationabsorbedby asurfaceisnecessary, whichneverthelessis notyetavailableintheliterature.Asan estimate,theresistanceforce againstdislocationabsorptionisproportional tothe rateofdislocationabsorption(measuredbytheplasticsliprateatthesurface).Moreover,theresistanceforceincreaseswith theaccumulateddensityofsurfacesteps.Thesecan beaccountedforinthe followingrate-dependentpower-lawrelation for

η

i^dis(α),

η

i^dis(α)=R^s_i^(α)

γ

˙^(α)

γ

˙0^s

ms

sgn

γ

˙^(α)

. (31)

γ

˙0^s isthereferencesliprateatthesurfaceandms isthesurfacerate-sensitivityexponent.ThesurfaceslipresistanceR^s_i^(α) isexpressedas

R^s_i^(α)=R^s₀^(α)+ β

H^s(αβ)

γ

acc^s(β), (32)

where

γ

acc^s(β)=

γ

˙^(β)

dt istheaccumulatedslipatthesurface,R^s₀^(α) denotestheinitialsurfaceslipresistance,andH^s(αβ) isthe dissipativesurfacehardening modulus. Forimperfect surfaceswith surfacecoatings,oxide layers orinitial defects, thesurfaceparametersdependontheinitialsurfacestate.Itisworthmentioningthatinsomegrainboundarymodels(van Beers etal., 2013), a similar formfor dissipativegrain boundary microforce is adopted, butthe slipresistance is simply assumedtobeconstant,bywhichthedissipativehardeningisignored.

Neglectingthemicrotraction_i^(α),themicroscopicboundaryconditionEq.(19)isrewrittenas

γ

˙^(α)=

γ

^˙0^s

₋

ξ

_i^(α)·Ni−

η

^en_i ^(α)

R^s_i^(α)

_m¹_s

sgn −

ξ

^(α)i ·N_i−

η

^eni ^(α)

(33)

whichconstitutesthegoverningequationfordislocation absorptionby surfaces.Bothenergeticanddissipativesurfaceef- fectsareaccountedfor.

3. Finiteelementimplementation

Forthe rate-dependentinitial boundary value problemfromthe presentrate-dependentmodel to be solved forarbi- trarygeometriesandboundaryconditions,adual-mixedfiniteelementformulation(Ekhetal.,2007)isadopted.Theforce balanceEq.(16) andtheevolutionEqs. (8)and(9)forthe GNDdensities arechosen asgoverningequations,and, hence, thedisplacement u andthe GNDdensities

ρ

^{ge(α) and}

ρ

^{gs(α) are considered as primary variables. The plasticslip}

γ

^{(α) is}

evaluated locallyat the integrationpoints. The treatment for thebalance of momentum (16)is straightforward, butdue tothesurfaceeffects,thedicretizationandimplementationoftheGNDevolutionequationsdeservesomeattention. Since theevolution equationsfor thetwo kinds ofdislocationsare similar, it isonly explainedfor edgeGND densities forthe sakeof brevity.Toobtain theweak formsofthegoverning equations,theyare multipliedby weighting functions

δ

^uand

δ ρ

^˙ge(α), whichvanishat thepartof thesurfacewithnaturalboundary conditions,andthen integratedoverthe volume.

Next,applyingthedivergencetheorem,weobtain

Fig. 1. Sketch of thin film benchmark problem.

Vi

[F^e·S^e]:

∇

i

δ

^udVi−

S^T_i

δ

^u·TidS^T_i =0, (34)

Vi

δ ρ

^˙_i^ge(α)

ρ

^˙_i^ge(α)−1

bs^(α)·

∇

i

δ ρ

^˙_i^ge(α)

γ

^˙(α)

dVi+

S^T_i

1

bs^(α)·Ni

δ ρ

^˙_i^ge(α)

γ

^˙(α)dSTi =0, (35) where the surface traction boundary condition(17) has been substituted into Eq.(34). The plastic slip rate

γ

˙^{(α) atthe} surfaceinEq.(35)isevaluated locallyvia themicroscopicsurfaceboundary condition(33),andtheinitial GNDdensities areneglected(butarestraightforwardtoincorporate).Forthespatialdiscretization,thevolumeofthebodyisdividedinto finite elementsineach ofwhich, followingthe standard Galerkinapproach,the unknown fieldsofthe displacement and GNDdensitiesandtheweightingfunctionsareapproximatedbytheirnodalvaluesmultipliedbyshapefunctions.Fortime integrationwithinatimeincrement[tn,t_n+1],theimplicitbackwardEulerschemeisadoptedsuchthat

γ

˙n^(α)+1=

γ

^(α)

^{t ,}

ρ

^˙n^ge+^(α)1 =

ρ

^ge(α)

^{t ,} F^p_n⁻₊₁¹=F^p_n⁻¹·

I−

α

γ

^(α)s(α)m^(α)

, (36)

where

γ

^(α)=

γ

n^(α)+1−

γ

n^(α),

ρ

^ge(α)=

ρ

^gen+^(α)1 −

ρ

n^ge(α), ^t=t_n+1−tn andthe subscript n identifies the quantities atthe nthtime step.Then, theweak forms(34)and(35)reduce toa highlynonlinear, stronglycoupledsystemofequations at eachtimestepwhicharesolvedbytheNewton–Raphsonmethod.

Fortheincorporationofthesurfacemodelintothefiniteelementformulation,theevaluationofsliprateatthesurface involvedinthesurfaceintegrationinEq.(35)isrequired.Tothisend,theincrementofplasticslip

γ

^{(α) atthesurfaceis}

determinedbythemicroscopicboundarycondition(33)as

γ

^(α)=

^t

γ

^˙0^s

₋

ξ

_i^(α)·Ni−

η

^en_i ^(α)

R^s_i^(α)

m¹s

sgn −

ξ

^(α)_i ·Ni−

η

^en_i ^(α)

, (37)

whereplasticslipsrequiredtoevaluatethetermsattherighthandsideareextrapolatedfromthebulk.Particularly,ifthe surfacedissipativemicroforceisignored,theincrementofplasticslip

γ

^{(α) atthesurfaceisevaluatedby}

−

ξ

i^(α)·N_i=

^(α)1

s^(α)·N_i

sgn

γ

^(α)+

γ

n^(α)

+

β

^(αβ)₂

s^(α)·N_i

m^(α)×N_i

s^(β)·N_i

m^(β)×N_i

γ

^(β)+

γ

n^(β)

, (38)

whichconstitutesthemicroscopicboundaryconditionsforapurelyenergeticsurface.

4. Numericalexample

Aplane strain benchmarkproblemofplaneconstrainedshearofa singlecrystallinethinfilmisstudied.As illustrated inFig.1,thefilmwiththicknesshinthex₂-directionisassumedtobeinfinitelylonginthex₁-direction,andwithoutloss ofgenerality, two active slipsystemswitheach ofthem definedby aslip directions^(α)=[cos

θ

α,sin

θ

α] anda slipplane normalm^(α)=[−sin

θ

α,cos

θ

α]areconsidered.Forthepresentplanestrainproblem,onlyedgeGNDsareinvolved.Thethin

film suffersa homogeneous external shear rate

γ

˙^{ext with its lower surfacebeing}constrained such that the macroscopic boundaryconditionsread

^u1

(

^x1,0

)

⁼

^u2

(

^x1,0

)

⁼

^u2

(

^x1,h

)

^=0,

^u1

(

^x1,h

)

⁼

γ

^˙exth

t, (39) where^u1 and^u2 denotethedisplacementincrementsaftertime increment^t.Microscopically,thesinglecrystalthin film isfree-standing, andthereforethe governing Eq.(33) fordislocation absorption by surfaces actsas the microscopic surfaceboundarycondition.Inaddition,tomodeltheinfiniteextensionofthefilminthex₁-direction,arepresentativepart ofthethinfilmwithlengthLisconsidered andperiodicboundaryconditionsare imposedforboththedisplacementand theGNDdensities

u1

(

^0,x2

)

^=u1

(

^L,x2

)

^{, u2}

(

^0,x2

)

^=u2

(

^L,x2

)

^,

ρ

^ge(1)

(

^0,x2

)

⁼

ρ

^ge(1)

(

^L,x2

)

^,

ρ

^ge(2)

(

^0,x2

)

⁼

ρ

^ge(2)

(

^L,x2

)

^{, (40)}

wheretheperiodiclength Lcan bearbitrarily chosen withoutaffectingthesolution.Forthespatialdiscretization,4-node bilinear quadrilateralelements with2×2full integrationare chosen forboth thedisplacements andthe GND densities.

Tonumericallyevaluatethesurfaceintegral,twoadditionalintegrationpointsareplacedonthesurfaceedges ofelements locatedatthesurface.

The elastic parameters and the magnitude of the Burgers vector representative of aluminium are taken here, i.e.,

μ

=26.3GPa,

ν

=0.33andb=0.286nm.Some ofthe plastic material parameters are takenasR^b₀^(α)=20MPa, H^b(αβ)= 200MPa,

γ

˙0^b=

γ

˙0^s=0.001s⁻¹,m_b=ms=0.05.The characteristic length scale isassumed to be constant, i.e., l=0.5

μ

^m.

The magnitudeofthe loading shearrate

γ

˙^{ext is 0}.001s⁻¹.To illustrate theinfluence ofenergetic anddissipativesurface effects,differentvalueswill be takenfortheremaining parameters involvedin thesurfacemodel,i.e., R^s₀^(α),H^s(αβ),^(α)₁ and₂^{(αβ)inthefollowing}simulations.Forsimplicity,thesurfaceparametersfromdifferentslipsystemsareassumedtobe thesame,i.e.,R^s₀^(α)=R^s₀,H^s(αβ)=H^s,1^(α)=1 and^(αβ)2 =2.

4.1. ValidationoftheFEMsolution

First,tovalidatethefiniteelement implementationsummarizedinSection3,theFEMsolutioniscomparedtothean- alyticalsolution of (Peng andHuang, 2017). The analytical solution in Pengand Huang(2017) only applies to the small deformationcasewithoutdissipativeeffects.However, ifthe appliedloadissmallanddissipativemicroforcesboth inthe bulk and the surfaceare ignored, the solution of the presentproblem should approach the analytical one. Thus, in the followingcomparison,thedissipativeforces

π

i^(α) and

η

^disi ^(α) aretemporarilyignored.Toavoidrepetition,theresultsfrom theanalyticalsolutionaredirectlygivenherewithoutpresentingthedetailedformulationwhichcanbefoundinPengand Huang(2017).Forboththenumericalandanalyticalcases,theassociatedparametersotherthanthosegivenabovearetaken ash=1

μ

^m,

θ

1=

π

/12,

θ

2=5

π

/12,1=10N/mand2=100N/m.InFig.2a,theshearstress

σ

12,whichishomogeneous alongthex₂-direction,isplottedasafunction oftheexternally appliedshear

γ

^{ext.Thenumericalresultscorrespondvery}

welltotheanalyticalresultsfor

γ

^ext<0.005.If

γ

^{extincreases,aslightdifferenceisnoticeablewhichisattributedtothefi-}

nitedeformationframeworkappliedhere.Thestress-straincurvesaredividedintothreestagesbytwoyieldingpointswhich denotetheonsetofdislocationabsorptionbysurfacesforthetwoslipsystems.Asdissipativeresistancesareignoredatthis stage,plasticdeformationoccursimmediatelyafterloading,and, hence,thereisnobulkyieldingpointonthestress-strain curve.InFig.2bandc,thedistributionsofplastic slipsandGNDdensities alongthex₂-direction areshown, wherethree valuesof

γ

^{extnamely0.001,0.002and0.005}representingthethreestagesrespectivelyaretaken.Forbothplasticslipsand GNDdensities,theFEMcurvescoincidealmostperfectlywiththeanalyticalcurves,indicatingtheexcellentmatchbetween theFEMandanalyticalsolutions.

4.2.Influenceofenergeticanddissipativesurfaceeffects

Inthepresentsurfacemodel,bothenergeticanddissipativesurfaceeffectsareinvolved.Here,toillustratetheirinfluence separately,two specific cases, i.e., the casewithonly energetic surfaceeffects and theone withonly dissipativesurface effectsareconsidered,respectively.Inbothcases,

θ

1=

π

/3and

θ

2=2

π

/3aretaken.

First,weconsiderthecasewithonlyenergeticsurfaceeffects.Fordifferentvaluesof1 and2,thestress-straincurve, theevolutionofthetotalaverageGNDdensity

ρ

¯^ge=h

0

ρ

^ge(1)

₊

ρ

^ge(2)

dx₂/h,theevolutionofsliprate

γ

˙^satthesur- face,andthedistributionofplasticslip

γ

^for

γ

^ext=0.01areplottedinFig.3a–d,respectively.Sincetheresultsforthetwo slipsystemsaresimilarduetosymmetry,onlythoseforslipsystem1aredisplayed.InFig.3a,thestress-straincurvesfor theenergeticsurfacemodelare generallydivided intothreestagesby two pointsrespectively denotingmacroscopicbulk yieldingandtheonsetofdislocationabsorptionbysurfaces(oralternatelysurfaceyielding).Similarstress-straincurveswith bothbulkandsurface/interfaceyieldingpointsarealsopredictedbysomeother modelswithsurface/interfaceeffects(see e.g.,Aifantisetal.,2006;vanBeersetal.,2013;FredrikssonandGudmundson,2005;Husseretal.,2017).Afterbulkyielding, thesurfaceinitiallyremainsimpenetrableduetothelackofdrivingforcefordislocationabsorption,andthehardeningrate isthesameasthatpredictedbymicrohardboundarycondition.Iftheexternallyappliedshearreachesacriticalvalue,the drivingforce issufficient to overcomethe energeticsurfaceresistance fordislocation absorption, surfaceyielding occurs,

Fig. 2. Comparison between FEM and analytical solutions: (a) stress-strain curve; (b) distribution of plastic slips; (c) distribution of GND densities. The FEM results correspond very well with the analytical results except the slight difference in stress-strain curve if γ^extis larger than 0.005 due to the finite-deformation framework in contrast to the small deformation theory employed in the analytical solution.

resultinginaconsiderabledecreaseinthehardeningrate.FromFig.3b,aftersurfaceyielding duetotheabsorptionofdis- locationsbythesurface,theaccumulationrateofthetotalaverageGNDdensitydecreasesinaccordancewiththedecrease ofhardeningrate.Comparingtheresultsfordifferentvaluesof1and2revealsthatthecriticalpointforsurfaceyielding onlydependsonthesurfaceenergyparameter1.Alarger1resultsinsurfaceyieldingatalarger

γ

^{ext.Thehardeningrate}

andtheaccumulationrateofGNDdensityaftersurfaceyieldingaregovernedbytheenergetichardeningmodulus2ina waythatalarger2 yieldsalarger hardeningrate.AsshowninFig.3c,aftersurfaceyielding,thesliprateatthesurface (measuringtherateofdislocationabsorptionbysurface)increasesimmediatelyfromzerotoaconstantleveldependingon thechosen magnitudeof2.Moreover,the surfaceresistanceagainst dislocationabsorption increaseswithincreasing 2. FromFig.3d,forsmaller1 and2,plasticslipsareeasiertodevelopanddistributemorehomogeneously.Thecurvespre- dictedbytheenergeticsurfacemodelareboundedbythoseformicrohardandmicrofreemodels,andtheenergeticsurface modelcanreadilyreduce tothemicrohard caseorthemicrofree caseby increasingthesurfaceparameters 1 and2 to infinityordecreasingthemtozero.

FromFig.3a,strainhardeningaftersurfaceyielding isinsignificantunlesstheenergetichardeningmodulusisimpracti- callylarge.²AsimilarconclusionwasdrawnfortheenergeticsurfacecontributioninthemodelofHurtadoandOrtiz(2012). There,animpracticallylargesurfaceparameterwasrequiredinorderforthemodeltocapturesize-dependentstrainhard- eningof micropillarsunder compression.Thecurrentresults implythat theenergetic surfaceeffectsare not sufficientto account forthe size-dependenthardeninginsmall-sizedsinglecrystals.Forthatreason, therole ofdissipativesurfaceef- fectsneedstobeinvestigated.Thus,weconsidernextthecasewithonlydissipativesurfaceeffects.Tothisend,theenergetic surfacemicroforcesareignoredbysetting1and2tozero.TheinitialsurfaceslipresistanceR^s₀andthedissipativesurface

2As estimated by Peng and Huang (2015) , the physically based energetic hardening modulus 2should be smaller than 100 N/m.

Fig. 3. Purely energetic surface model: (a) stress-strain curves; (b) evolution of the total average GND density; (c) evolution of slip rate at the surface; (d) distribution of plastic slip for γ^ext= 0 . 01 . Surface yielding depends on 1. Surface hardening, the accumulation rate of the total average GND density after surface yielding, and the slip rate at the surface at the steady-state are governed by the energetic surface hardening modulus 2. For smaller 1and 2, plastic slips develop faster and distribute more homogeneously after surface yielding. The results of the energetic surface model are bounded by those of microhard and microfree models.

hardeningmodulusH^sarevaried. InFig.4a–d,thestress-straincurve,theevolutionofthetotalaverageGNDdensity

ρ

¯^ge, theevolutionofsliprate

γ

˙^{s atthesurface,andthe}distributionofplasticslip

γ

^for

γ

^ext=0.01 fordifferentvaluesofR^s₀ andH^saredisplayed,respectively.Thepurelydissipativesurfacemodelpredictsplasticbehaviorquitesimilartothatofthe purelyenergetic one.Particularly, thesurfaceyielding dependsonthe initialslipresistance R^s₀,while thehardeningrate, theaccumulation rateofthetotal GNDdensityafter surfaceyielding, andthesliprateat thesurfaceatsteady-state are governedbythedissipativesurfacehardeningmodulusH^s.

Inviewofreversibilityassociatedwithsurfacestepformation,weexploretheinfluenceofsurfaceeffectsontheplastic behavior duringunloading. Tothisend, the stress-straincurve, the evolution ofsliprate atthe surfaceduringa loading andunloading cycle predictedby the purely energetic surface model,the purely dissipative surfacemodel andthe gen- eralsurfacemodel withboth energetic anddissipativesurface effectsare displayed inFig. 5a andb,respectively, where 1=10N/m,2=100N/m,R^s₀=5N/mandH^s=200N/m.The loading isapplied up to

γ

^ext=0.01,followed by full un- loading.Fromthefigures,thecurvespredictedbythetwomodelsaresimilarduringloading.Forunloading,thesituationis different.Forthepurelyenergeticsurfacemodel,duringunloading,theenergeticsurfacemicroforce,whichcanberegarded asasurfaceback-stress,actsasthedrivingforcefortherecovery ofsurfacesteps(inversetodislocationabsorptionduring loading).Hence,bulkyieldingandsurfaceyieldingoccuratthesametime.Thehardeningrateandthesliprateatthesur- faceimmediatelyafteryielding (thesecondunloadingstage)areidenticaltothatofthethirdloadingstagewithbothbulk andsurfacehardening.If

γ

^{extdecreasestoacriticalvalue,theenergetic}microforcesdecreasetozero,thesurfacebecomes impenetrableagain. After that the hardening rate remains the same asthat ofthe second loading stage with onlybulk hardening,andthe sliprateat thesurfacedecreases to zero.Forthe purely dissipativesurfacemodel,during unloading, duetotheexistenceofthedissipativesurfacemicroforce,therecoveryofsurfacestepsdoesnottakeplaceuntilthedriving

Fig. 4. Purely dissipative surface model: (a) stress-strain curves; (b) evolution of the total average GND density; (c) evolution of slip rate at the surface;

(d) distribution of plastic slip for γ^{ext= 0 . 01} . Surface yielding depends on initial surface resistance R ^s₀. Surface hardening, the accumulation rate of the total average GND density after surface yielding, and the slip rate at the surface at the steady-state are governed by dissipative surface hardening modulus H ^s. For smaller R ^s₀and H ^s, plastic slips develop faster and distribute more homogeneously after surface yielding. The results of the dissipative surface are bounded by those of the microhard and microfree models.

force (the net contributionofbulk andsurfaceenergetic microforces)is sufficient toovercomethe dissipativeresistance.

Hence,thebulkyieldsfirstand,subsequently,thesurfaceyieldingoccursonlyifthedrivingforceovercomestheresistance.

Thehardeningrateandthesliprateatthesurfaceforthesecondandthirdunloadingstagesareidenticaltothoseforthe correspondingloadingstages.Forthegeneralsurfacemodelwithbothdissipativeandenergeticsurfaceeffects,theloading andunloadingbehaviorsaresimilartothoseofthepurelydissipativesurfacemodel.Inaddition,apronouncedBauschinger effectisalsoobservedinFig.5a.

4.3. Influencesoffilmthicknessandorientationofslipsystems

Toexploretheinfluencesofthefilmthicknessh,fordifferentvaluesofh,thestress-straincurve,theevolutionofthetotal averageGNDdensity

ρ

^{ge ,theevolutionofsliprate}

γ

^{˙ atthesurface, andthe}distributionofplasticslip

γ

^for

γ

^ext=0.01 are plottedin Fig.6a–d, where

θ

1=

π

/3 and

θ

2=2

π

/3,1=5N/m,2=100N/m,R^s₀=5N/mand H^s=200N/m. The resultsdisplayedinFig.6clearlydepict asize-dependentmaterialresponse.Inthinnerfilms,theflowstress andthetotal averageGNDdensityataspecific

γ

^{extare larger,whichindicatesthat “smallerisstronger”.AsshowninFig.6aandc,in}

thinnerfilms, surfaceyieldingoccursatasmallerapplied shear

γ

^{ext,andthesliprateatthesurface(measuringtherate}

of dislocation absorption)atthe steady state islarger, which impliesthat surface yielding occursmore easily atsmaller scales.Infact,intheexperimentsby Gruberetal.(2008),itisfoundthat forsingle-crystallineAuthinfilms withsmaller thickness,afterbulkyielding,thestrainhardeningratesdecreasesuddenlyataspecificpointwhichissuspectedtobethe surfaceyieldingpoint,whileforfilms withlargerthickness, nosuchphenomenonoccurs, whichisqualitativelyconsistent withthepredictedresultshere.

Fig. 5. Purely dissipative surface model versus purely energetic surface model: (a) stress-strain curves; (b) evolution of slip rate at the surface. A pro- nounced Bauschinger effect of earlier re-yielding is observed. The results predicted by the purely energetic and dissipative surface models are similar during loading but differ during unloading.

Fig. 6. Size effect: (a) stress-strain curves; (b) evolution of the total average GND density; (c) evolution of slip rate at the surface; (d) distribution of plastic slip for γ^ext= 0 . 01 . In thinner films, surface yielding occurs earlier. A smaller thickness results in a larger flow stress and a larger total average GND density at a specific load, and a larger surface slip rate at the steady-state. In the middle of the film, a smaller thickness yields a smaller plastic slip. Near the surface, a smaller thickness gives a larger plastic slip.