Multiple slip dislocation patterning in a dislocation

0 downloads 0 Views 5MB Size Report
Oct 10, 2017 - eff is called effective shear stress, tp is the Peierls stress and the threshold stress ta th is given by dislocation interactions on various slip ...
International Journal of Plasticity 100 (2018) 104e121

Contents lists available at ScienceDirect

International Journal of Plasticity journal homepage: www.elsevier.com/locate/ijplas

Multiple slip dislocation patterning in a dislocation-based crystal plasticity finite element method € bes c, d, D. Raabe c N. Grilli a, b, *, K.G.F. Janssens a, J. Nellessen c, S. Sandlo a

Laboratory for Nuclear Materials, Nuclear Energy and Safety Department, Paul Scherrer Institut, 5232 Villigen PSI, Switzerland  NXMM Laboratory, IMX, Ecole Polytechnique F ed erale de Lausanne, 1015 Lausanne, Switzerland c Max-Planck-Institut für Eisenforschung GmbH, Department for Microstructure Physics and Alloy Design, 40237 Düsseldorf, Germany d Institute of Physical Metallurgy and Metal Physics, RWTH, Aachen University, D-52056 Aachen, Germany b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 4 May 2016 Received in revised form 15 September 2017 Accepted 29 September 2017 Available online 10 October 2017

Dislocation structures forming during cyclic loading of fcc metals are fatigue damage precursors. Their specific structures are caused by the motion and interactions of dislocations. Depending on the load conditions, the grain orientation, the stacking fault energy, a variety of different dislocation structures appear in the material such as labyrinths, cells, veins and persistent slip bands. We present a continuum dislocation-based model for cyclic fatigue and incorporate it into a crystal plasticity finite element solver. A method for the simulation of dislocation junction formation is introduced, which reproduces the behaviour of discrete objects, such as dislocations, in a continuum framework. The formation of dislocation walls after 50 and 100 deformation cycles at 0.95% and 0.65% strain amplitude starting from an initial random dislocation distribution is predicted for 〈001〉 and 〈110〉 oriented crystals. Simulations and cyclic tension-compression experiments of polycrystalline 316L stainless steel are performed to compare our model with another model based on edge and screw dislocation densities. The simulated dislocation structures and experimental results, obtained with the electron channeling contrast imaging technique, are compared using a 2D orientation distribution function of the dislocation structures. The dominant orientation of dislocation walls is predicted by the new model; it turns out to be perpendicular to the intersection line between the two slip planes involved in their formation and at an angle of around 45o from the loading axis. This agrees well with the experimental observations and represents a step forward for understanding the formation mechanism of these dislocation structures. © 2017 Elsevier Ltd. All rights reserved.

Keywords: Dislocations Dislocation-based crystal plasticity Fatigue Finite element method

1. Introduction Cyclic fatigue in fcc metals is characterised by the formation of dislocation structures, which are localized regions where the dislocation density becomes higher than in pristine materials (Li et al., 2011). Dislocation structures are constituted of dislocation walls separating low dislocation density regions known as channels. The appearance of these dislocation structures is associated with hardening of metals and precedes crack formation (Pham and Holdsworth, 2012). Phenomenological models are often used to describe damage during the first stage of fatigue (Pirondi et al., 2006); however, dislocation-based

* Corresponding author. Laboratory for Nuclear Materials, Nuclear Energy and Safety Department, Paul Scherrer Institut, 5232 Villigen PSI, Switzerland. E-mail address: [email protected] (N. Grilli). https://doi.org/10.1016/j.ijplas.2017.09.015 0749-6419/© 2017 Elsevier Ltd. All rights reserved.

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

105

models at the sub-micrometer length scale contain information related to physical parameters and the microstructure, hence, they are helpful in enlightening the underlying mechanisms (Korsunsky et al., 2012). Current research on models for dislocation structures formation is mostly focused on monotonic loads Sandfeld and Zaiser, (2015), for which the dislocation cell structure has been found (Xia and El-Azab, 2015). Cyclic simulations of single crystals oriented for single slip, for which vein-channel structures and persistent slip bands (PSBs) are observed Li et al. (2011), have pre s et al. (2008); Fivel (2008) computational been carried out using both continuum Pontes et al. (2006) and discrete De approaches starting from random initial dislocation distributions. The use of discrete dislocation dynamics simulations for fatigue is limited by the computational resources needed to reach few cycles, and is restricted to simplified boundary conditions (Kubin, 1002). Details about the positions of single dislocations are not a necessary condition to study the micrometerscale features of dislocation structure formation, therefore continuum models, based on dislocation densities, are a suitable alternative. A finite element method based approach is capable to take into account displacement controlled boundary conditions and the internal stresses caused by plastic deformation. For this reason, we use the crystal plasticity finite element (CPFE) method Roters et al. (2012), a computational method frequently used to simulate polycrystals, mostly considering monotonic loads only (Leung et al., 2015). The CPFE method needs to be complemented by constitutive equations to calculate the plastic strain rate of every slip system. These equations can be based on a phenomenological (Manonukul and Dunne, 2004) or, as in the present work, on a dislocation dynamics model. In polycrystalline samples the specific arrangement of grains and their orientation determines the strain and stress states of a material at the micrometer scale, which controls the active slip systems according to Schmid's law (Roters et al., 2010). The CPFE method has been proven suitable for small-scale plasticity simulations, because the plastic deformation of all the slip systems is considered and, therefore, the active ones can be predicted (Maab et al., 2009). The interaction of dislocations on different active slip systems determines the orientation of dislocation walls. Monotonic deformation experiments show 3D cell structures at 2% strain amplitude (Jakobsen et al., 2006). The activation of cross slip is associated with the appearance of these structures, which do not emerge if this mechanism is suppressed in 3D discrete dislocation dynamics simulations (Madec et al., 2002). The dislocation walls during monotonic deformation tend to align along the intersection of cross-slip planes (Xia and El-Azab, 2015). Cyclic deformation experiments show a variety of differently oriented dislocation structures in fcc metals, such as copper, nickel and silver (Li et al., 2011): PSBs and vein-channel structures coexist in 〈011〉 oriented crystals and walls perpendicular to the Burgers vector of the most active slip system have also been observed (Li et al., 2011). Cell structures, composed mostly of coplanar dislocations, were reported to form in 〈111〉 oriented crystals Li et al. (2009a), and labyrinth structures in 〈001〉 oriented crystals. The wall orientation of labyrinth structures in ½001 oriented crystals is typically perpendicular to ½001 and ½100 axes Li et al. (2011), thus perpendicular to the sum of the Burgers vectors of the primary and the secondary critical slip systems. In 〈001〉 oriented crystals, when the conjugate and the critical slip systems are activated, Lomer-Cottrell and Hirth junctions are formed, respectively. The Lomer-Cottrell junction is known to be stronger than the Hirth junctions. This means that, in presence of Lomer-Cottrell junctions, a higher stress is required to activate the subsequent dislocation movement. Therefore, the dislocation motion and multiplication is slowed down in regions where the primary and conjugate slip systems are activated (Jin and Winter 1984a). This is not the case for the Hirth junction, which does not retard the contemporary motion of dislocations belonging to the primary and to the critical slip system (Jin, 1983; Jin and Winter 1984b; Li et al., 2009b). This suggests that the activation of the critical slip system, and consequently the formation of Hirth locks, is important for the labyrinth structure formation (Li et al., 2011). Further experimental investigations have found wall orientations intermediate between f100g and f210g L'Espe'rance et al. (1986), explained with low energy arrangements of dislocation loops (Dickson et al., 1986). The collinear reaction, which is the interaction of a slip system with its cross-slip system, is also important to determine the geometrical configuration of dislocations during loading (Madec et al., 2003). However, active slip systems interacting by the collinear interaction are not commonly observed in copper (Devincre et al., 2005, 2006). There have been only a few attempts to model dislocation structure formation under multiple slip conditions. Approaches based on reaction-diffusion equations Aoyagi et al. (2013) have been applied to simulate the formation of labyrinth structures Pontes et al. (2006), where mobile and immobile, positive and negative edge dislocations are the state variables. These attempts successfully reproduced two-dimensional dislocation patterning, resulting in dislocation walls perpendicular to the Burgers vectors of the two active slip systems. This is different from experimental observations, where these walls are oriented perpendicular to the sum and difference of both Burgers vectors. This discrepancy is due to the flux terms used to model the motion of edge dislocation segments along the two slip directions. Opposite signed edge dislocations, belonging to the two active slip systems, move away from emerging channels along the Burgers vector direction, interact and form dislocation dipoles. Therefore, dislocation walls do not align with the dislocation junction direction and remain oriented perpendicular to the Burgers vectors. Thus, it appears necessary to introduce dislocation junctions and their specific geometrical configuration in constitutive models, as has been done in discrete dislocation dynamics models Martínez et al. (2008) but not in continuum formulations. Testing these constitutive models on 〈001〉 oriented crystals can indicate if the dislocation junction mechanism is sufficient to describe labyrinth structures. In this paper we apply the crystal plasticity finite element (CPFE) method and dislocation-based constitutive laws, implemented as a DAMASK subroutine Roters et al. (2012), to model dislocation structures in polycrystalline fcc metals. The main novelty introduced is the treatment of dislocation junction formation in a continuum framework, based on the linear combination of edge and screw dislocation densities into a density aligned with dislocation junctions. When applied to the simulation of cyclic plasticity in 〈001〉 oriented crystals, the new model is capable to predict the correct orientation of

106

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

dislocation walls, shedding light on the formation mechanism of these dislocation structures. The simulations are analysed using a 2D orientation distribution function Gasparyan and Ohanyan (2015) and compared to experimentally imaged dislocation structures after cyclic deformation using the electron channeling contrast imaging (ECCI) technique Zaefferer and Elhami (2014) on 316L austenitic steel specimens. We show that previous theories Dickson et al. (1986), based on low energy dislocation structures, are incomplete and need further modifications to correctly predict the orientations of dislocation structures. This paper is organised as follows: in section 2 the CPFE framework is presented; in section 3.1 and 3.2 two constitutive models, one including dislocation junctions and another one based on edge and screw dislocations, are introduced; in section 4.1 a proof of concept is shown, where two dislocation loops interact and form a single dislocation junction; in section 4.2 polycrystal simulations of three different grains using the two constitutive models are shown; in section 4.3 the simulated hardening curves are reported and compared with the experimental results; in section 5 simulation results are discussed and conclusions concerning the dynamics of dislocation structures formation are drawn. 2. CPFE method and dislocation-based constitutive models The crystal plasticity finite element method combines a finite element solver with a micro-scale description of plasticity in the material. Mesh and boundary conditions are introduced in a finite element model. The deformation gradient at every element is decomposed into an elastic and plastic part (Roters et al., 2012):

F ¼ Fe ,Fp :

(1)

The time evolution of the plastic part Fp is determined by the motion of dislocations, which belong to the 12 slip systems of fcc materials, according to Orowan's law:

F_ p ,F1 p ¼

12 X a ba ; m 5n g_ ap c

(2)

a¼1

g_ ap ¼ ra va b ;

(3)

where b is the Burgers vector, va and ra are the dislocation velocity and total dislocation density on the slip system a. The a b a are the slip direction and normal of the slip system a. The stress tensor (second Piola-Kirchhoff stress) is vectors c m and n calculated from the elastic part of the deformation gradient:



 ℂ T Fe ,Fe  I ; 2

(4)

where ℂ is the stiffness tensor. The stress tensor respects the equilibrium conditions and a Newton-Raphson scheme is used to calculate the plastic strain gradient (Roters et al., 2012). This method is able to capture the elastic field created by dislocations and their mutual elastic interactions. ! The time evolution of the dislocation density state vector rð r ; tÞ, whose components represent densities of different dislocation types, is given by rate equations:

r_ ¼ gðrð! r ; tÞÞ ;

(5)

where the function g describes the dislocation processes in the material (Roters et al., 2010). In the following edge and screw dislocation densities will be indicated by the symbols r ⊥ and r N , the symbol r ⊥N is used in equations that are valid for both edge and screw dislocations. Hereafter the dislocation processes common for both models are described, namely annihilation and flux. The annihilation rate is proportional to the meeting frequency of positive and negative signed dislocations (edges and screws) (Ma and Roters, 2004). For instance the annihilation rate of positive and negative dislocation densities on a slip system a is then given by:





 v ⊥N  ; r_ a⊥Nþ;ann ¼ r_ a⊥N;ann ¼ 4 bd ⊥N ra⊥Nþ ra⊥N ! a

(6) a

where b d ⊥N is the annihilation distance of edge and screw dislocations and ! v ⊥N is the dislocation velocity on the slip system a. The motion of dislocations between neighbouring elements is given by a flux term:

 a  J a⊥N ¼ V, ra⊥N ! v ⊥N :

(7)

a Here the velocity ! v ⊥N is directed perpendicular to the line direction of edge or screw dislocations, its magnitude va is assumed isotropic and proportional to the resolved shear stress on each slip system (Sandfeld and Zaiser, 2015):

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

107

   va ¼ B ta   tath  tp ¼ Btaeff ;

(8)

where B is the drag coefficient, taeff is called effective shear stress, tp is the Peierls stress and the threshold stress tath is given by dislocation interactions on various slip systems (Franciosi and Zaoui, 1982):

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 12 uX tath ¼ mbt Gab rb ;

(9)

b¼1

where m is the shear modulus and Gab is the matrix representing the strength of interactions between dislocations of slip systems a and b (Arsenlis and Parks, 2002). The dislocation velocity va in (8) is non-zero only if the effective shear stress taeff is positive. The sign of the dislocation velocity is the same as the one of ta . A temporary decrease of the total dislocation density in certain elements during the formation of dislocation structures is observed because of dislocation fluxes, leading to very high values of the dislocation velocity, as determined by (3). For this reason a maximum velocity vmax is set in the model to avoid numerical unconvergences without decreasing the computational time step. The typical effective resolved shear stress reached in our simulations is around 20 MPa, therefore a value vmax ¼ 1 mm=s has been chosen, which is consistent with the measured velocity in single crystals for that stress value (Turner and Vreeland, 1970). The time step in the simulations is dt ¼ 2,105 s to keep the travelling distance of a dislocation in one time step smaller than the element size. Model parameters mentioned above are included in Table 1. 3. Dislocation-based constitutive models 3.1. Junction constitutive model The new model uses multiple dislocation densities to differentiate between the orientations of a dislocation line. These state variables can be defined starting from the higher order dislocation density ra ðp; 4Þ of the continuum dislocation dynamics theory Hochrainer et al. (2014), which quantifies the density of dislocations at a point p on the slip system a with line direction oriented at an angle 4 with respect to the Burgers vector. Integration of ra ðp; 4Þ over prescribed angular intervals as defined in the following, leads to different dislocation densities. In case two slip systems are active, for which the slip plane ! ! b1  n b 2 Þ, as shown in Fig. 1 (a). normals are n1 and n2 , the dislocation junction vector l lock can be defined as l lock ¼ ð n ! Dislocation segments orthogonal to l lock are included in the dislocation densities ra£þ and ra£ , as shown in Fig. 1 (a). This orientation corresponds to dislocations gliding inside channels between neighbouring dislocation walls formed by Hirth locks in labyrinth dislocation structures (Li et al., 2011). The positive or negative sign depends on the position of the dislocation line in a clockwise oriented dislocation loop. Similarly ra==þ and ra== are defined as the densities of dislocation segments parallel to ! l lock . These dislocations have the same orientation as dislocation junctions. General dislocations, with an orientation intermediate between £ and ==, are included in other dislocation densities. For instance by integrating ra ðp; 4Þ over the angular interval Dq=þ;£ shown in Fig. 1 (a), which includes orientations intermediate between positive == and negative £, the density ra==þ;£ can be defined. The 8 different densities used in the model are listed in Table 2. If more than two slip systems are ! active, a fixed reference slip system is chosen to define the dislocation junction vector l lock, r== and r£ dislocations. The velocity vectors of different dislocation densities are orthogonal to the corresponding dislocation lines; for instance 1 1 1 the velocity vectors ! v £þ , ! v £ and ! v ==þ are depicted in Fig. 1 (b). Other dislocation densities have velocity vectors oriented at ! a a 45 or 135 , respectively, with l lock ; for instance ra==þ;£þ has velocity vector ð! v ==þ þ ! v £þ Þ. This definition of dislocation densities is suitable for modelling the dynamics of dislocation junctions in labyrinth ! dislocation structures, since under an applied stress the dislocations ra£þ and ra£ move along l lock , as shown in Fig. 1 (b), a a while the dislocations r==þ and r== can interact and form immobile junctions. A dislocation loop has to remain connected, therefore a new dislocation line is created that lengthens the dislocation junction. This dislocation multiplication process is accounted for by the dislocation curvature density Hochrainer et al. (2007), which concentrates between == and £ dislocation segments, for instance in r1==þ;£þ in Fig. 1 (a). As a consequence the densities ra==þ , ra== , ra£þ and ra£ do not participate in the dislocation multiplication process because they are constituted of dislocation segments with zero curvature. A curved dislocation segment, such as ra==þ;£þ , is always connected  toaone  == and one £ segment and its motion creates new == a  dislocation lines at a rate given by the velocities ! v £þ  and ! v £ . The number of these curved segments per unit volume is Table 1 Model parameters (Turner and Vreeland, 1970; Pham et al., 2013; Catalao et al., 2005; Sauzay, 2009; Misra et al., 2010) and elastic constants (Reed and €chen et al., 2008) of 316L steel. Horiuchi, 1983; Jo b

b d⊥

b dN

B

tp

vmax

0.258 nm

1.3 nm

50 nm

100 (mm/s)/MPa

0.112 MPa

1000 mm/s

C11

C12

C44

m

206 GPa

133 GPa

119 GPa

86 GPa

108

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

! 1 1 1 Fig. 1. (a) Definition of the dislocation junction vector l lock and of the different dislocation densities. (b) Velocity vectors ! v £þ , ! v £ and ! v ==þ .

Table 2 Dislocation densities in the junction constitutive model.

ra==þ ra==þ;£þ

ra== ra==þ;£

ra£þ ra==;£þ

ra£ ra==;£

given by the ratio between the corresponding density and a fixed average segment length L, for instance ra==þ;£þ ==L. The dislocation multiplication law can be written as:



_a

r==þ;mul ¼



 v £þ  ra==þ;£þ ! a

L



þ



 v £  ra==þ;£ ! L

a

:

(10)

Similar equations for r_ a==;mul, r_ a£þ;mul and r_ a£;mul can be obtained by substituting the curved densities in (10) ra==þ;£þ and ra==þ;£ with the two curved densities with ==, £þ or £ into the subscript. This dislocation multiplication law satisfies the

connectivity of dislocations Leung and Ngan (2016) because it is obtained by integrating over prescribed angular intervals the higher order dislocation density ra ðp; 4Þ Hochrainer et al. (2014), which respects connectivity; however, it can model only convex segments with a constant average curvature L. The value of L is 0.1 mm, so that a curved dislocation segment has a length corresponding to the distance between two neighbouring integration points. Specific equations for cross slip are introduced to model the motion of screw dislocations on secondary slip planes and their consequent multiplication as Frank-Read sources (Bitzek et al., 2008; Messerschmidt and Bartsch, 2003), as shown in Fig. 2 (a). This process generates new curved dislocation segments at a rate proportional to the cross slip rate:

 t Vact ; R ¼ nexp  III kB T

(11)

where n ¼ 2,1015 s1 is the attack frequency, determined from molecular dynamics simulations Vegge et al. (2000), tIII is the pre s et al., 2008; De pre s, 2004). The stage three stress and Vact is the activation volume. The value of tIII is around 56 MPa (De activation volume has been measured to be around 400 atomic volumes (Bonneville et al., 1988). The atomic volume of 316L pre s et al. (2008) and, therefore, Vact ¼ 5,1027 m3. Room temperature is used in all following stainless steel has been used De simulations. During every cross slip event a screw dislocation segment of length wz44b Bonneville et al. (1988) transfers on a parallel slip plane. The curved segment expands, as shown in Fig. 2 and a dislocation loop, whose average length is L, forms. Therefore, the number of loops created per unit volume during a time interval dt is:

  raNþ þ raN Rdt dracs ; ¼ w w

(12)

where dracs is the density of screw dislocation segments on the slip system a that cross-slipped. The rate at which new curved dislocation density is created is:

r_ a==þ;£þ;crossslip ¼

r_ acs L w

¼

 RL  a r þ raN ; w Nþ

and similar equations hold for the other curved densities.

(13)

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

109

Fig. 2. (a) The process of cross slip of a screw dislocation segment from one slip plane to another. (b) The resulting Frank-Read source generates new curved dislocation segments.

The screw dislocation densities raNþ and raN are obtained by projecting == and £ dislocation segments along the Burgers vector:

raNþ ¼ ra£þ cosF þ ra==þ sinF ;

(14)

raN ¼ ra£ cosF þ ra== sinF ;

(15)

where 0  F  p=2 is the angle represented in Fig. 1 (b):

 ! !   l  lock,b  ; sin F ¼  !  !  l lock , b 

(16)

! !  ! !     where  l lock  and  b  are the norms of l lock and b . The annihilation rate is proportional to the product of opposite signed edge and screw dislocation densities, obtained by projecting == and £ dislocation segments. For instance, the annihilation rate of r==þ dislocations is given by:









    a   a  v N  4b v N : d ⊥ ra⊥ , ra=þ cosF ,! r_ a==þ;ann;screw ¼ 4 bd N raN , ra=þ sinF ,!

(17)

The two terms in (17) take into account the annihilation of screw and edge parts of r==þ dislocation segments. 3.2. Edge-screw constitutive model In this section the edge-screw model Arsenlis and Parks (2002), used for comparison with the junction model, is introduced. Positive and negative, edge and screw dislocation densities (r⊥þ , r⊥ , rNþ and rN ) can be defined starting from the higher order dislocation density ra ðp; 4Þ. This is the density of dislocations at a point p with the line direction at an angle 4 to the Burgers vector, considering segments orthogonal and parallel to the Burgers vector (Arsenlis and Parks, 2002). Dislocations with an orientation intermediate between edge and screw are included in other dislocation densities; for instance r⊥þ;Nþ is defined by integrating ra ðp; 4Þ over orientation angles 4 between the positive edge and the positive screw dislocations, in analogy with the dislocation junction model. The 8 different densities used in the model are listed in Table 3. This definition of dislocation densities is suitable for modelling single slip dislocation structures, where edge dislocation walls are separated by channels (Laird et al., 1986). In this configuration, screw dislocations glide in channels while edge dislocations interact and form immobile dipoles. Thus the curvature density is concentrated between edge and screw dislocation segments (Grilli et al., 2015), in analogy with the geometry shown in Fig. 1 (b) for the dislocation junction model. The dislocation multiplication law for positive edge dislocations can be written as:

Table 3 Dislocation densities in the edge-screw model.

ra⊥þ ra⊥þ;Nþ

ra⊥ ra⊥þ;N

raNþ ra⊥;Nþ

raN ra⊥;N

110

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121 Table 4 Velocity vector directions of the different dislocation densities, referred to the coordinates of the slip system a in Fig. 3. a ! v ==þ ¼ va ðcosF; sinFÞ a ! v £þ ¼ va ðsinF; cosFÞ      a ! v ==þ;£þ ¼ va cos F þ p4 ; sin F þ p4      a ! v ==;£þ ¼ va  cos F  p4 ; sin F  p4 a ! v ⊥þ ¼ va ð1; 0Þ a ! v Nþ ¼ va ð0; 1Þ pffiffiffi pffiffiffi a ! v ⊥þ;Nþ ¼ va ð1= p 2 ;ffiffiffi1= p 2 Þffiffiffi a ! a v ¼ v ð1= 2 ; 1= 2 Þ

a

! v == ¼ va ðcosF; sinFÞ a ! v £ ¼ va ðsin Þ  F; cosF   a ! ¼ va cos F  p ; sin F  p v ==þ;£

4

4

     a ! v ==;£ ¼ va  cos F þ p4 ; sin F þ p4

a ! v ⊥ ¼ va ð1; 0Þ a ! v N ¼ va ð0; 1Þpffiffiffi pffiffiffi a ! v ⊥þ;N ¼ va ð1= p 2 ;ffiffiffi1= p 2 Þffiffiffi a ! a v ¼ v ð1= 2 ; 1= 2 Þ

⊥;Nþ

⊥;N

! ! Fig. 3. Coordinate system on the slip plane a to define the velocity vector directions: l lock is the dislocation junction vector and b the Burgers vector.



_a

r⊥þ;mul ¼



 v Nþ  ra⊥þ;Nþ ! a

L



þ



 v N  ra⊥þ;N ! L

a

;

(18)

and similar equations hold for r_ a⊥;mul, r_ aNþ;mul and r_ aN;mul . Cross slip is modelled as in section 3.1 and the cross slip rate for r⊥þ;Nþ can be written as in (13):

r_ a⊥þ;Nþ;crossslip ¼

r_ acs L w

¼

 RL  a rNþ þ raN : w

(19)

Projection equations (14) and (15) are not needed in the edge-screw model because positive and negative screw dislocation densities raNþ and raN are state variables. The velocity vectors of the different dislocation densities on a slip system a are summarized in Table 4. 4. Results 4.1. Two dislocation loops forming a dislocation junction To illustrate the junction model in the presence of dislocations belonging to two different slip systems, we simulate the motion of two dislocation loops on intersecting slip planes. The model geometry is constituted of two planes with a thickness of 200 nm, corresponding to one single element, that intersect along a line of elements, as shown in Fig. 4 (a). These planes correspond to the ð111Þ and ð111Þ slip planes of the FCC crystal and the angle between them is around 109o . A displacement along the z axis is imposed on the upper surfaces of the two planes while the lower surfaces are held fixed; the resulting z component of the displacement is shown in Fig. 4 (a). Using these boundary conditions, the resolved shear stress on the two slip systems is uniform. Fig. 4 (b) shows the evolution of the total dislocation density r on the two slip systems during changes in the resolved shear strain g. The initial dislocation density of the two loops is set at 1014 m2, thus each integration point, having a distance of 100 nm, contains the equivalent of a single dislocation. The two dislocation loops expand until they meet at the intersection of the two planes, where a dislocation junction is formed. The ra==þ dislocation density becomes immobile due to the dislocation interaction equation (9); without this interaction the dislocation density would exit the geometry because of the open boundary conditions for the dislocation fluxes, as shown in Fig. 4 (c). The sections of the loops that do not form the junction expand further and a new dislocation line is created at the connections between == and £ segments, as predicted by the multiplication law (10). The flux of the dislocation density from one element to its neighbour and the elastic interactions, leading to repulsion of same-sign dislocations Leung and Ngan (2016), cause the spread in width of the dislocation density, as shown in Fig. 4 (b). The stacking-fault between partial dislocations should be modelled to limit this core spreading, but the element size used in the present simulations (200 nm) is not suitable to resolve it. Another strategy to limit the core spreading has been recently proposed Leung and Ngan (2016), which is based on the relocation to nearby cells of the dislocation density below a certain threshold. This simulation is made without the cross slip process described by (13).

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

111

4.2. Dislocation patterning Both models are applied to the simulation of fatigue experiments on polycrystalline 316 austenitic stainless steel, where the electron channeling contrast imaging (ECCI) technique Nellessen et al. (2015) is employed. Bone-shaped specimens are characterized using EBSD and then deformed cyclically up to 100 cycles using a tensile-compression module under displacement control with a strain rate of 0:5,103 s1. In situ DIC measurements are used to determine the local strain distribution. The measured strain amplitude is averaged over regions of 500 mm, the same order of magnitude of the simulation representative volume size. Three grains, whose orientations are close to ½3 29 1, ½1 19 0 and ½9 10 0, are selected and analyzed by ECCI, which enables the observation of dislocations with different Burgers vectors. Evidence of oriented dislocation walls are provided by electron channeling contrast imaging observations, as shown in Figs. 5 (a), 6 (a) and 7 (a). EBSD data are used to build the corresponding polycrystal geometries in Figs. 5 (b), 6 (b) and 7 (b), whose dimensions are 186 mm  203 mm  6 mm, 182 mm  195 mm  6 mm and 305 mm  306 mm  6 mm, and to assign the proper stiffness tensor to every grain depending on the orientation. The mesh used (Figs. 5 (b), 6 (b) and 7 (b)) is coarser in the grains surrounding the central one and refined in the centre, where 200 nm elements are generated inside a parallelepiped with dimensions 4 mm  4 mm and depth 0.6 mm, close to the upper surface. The strain in the central grain is determined by the polycrystalline

Fig. 4. (a) Displacement field along the z axis. (b) Evolution of the total dislocation density r of two dislocation loops forming a dislocation junction at the intersection of their slip planes. (c) Total dislocation density after deformation without the interaction term ðtath ¼ 0Þ.

Fig. 5. (a) ECCI image taken in the centre of the ½3 29 1-oriented grain after 50 cycles; ð111Þ and ð111Þ slip plane traces are visible. (b) Simulation set-up of the polycrystalline structure (LD: loading direction).

112

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

Fig. 6. (a) ECCI image taken in the centre of the ½1 19 0-oriented grain after 100 cycles; ð111Þ, ð111Þ and ð111Þ slip plane traces are visible. (b) Simulation set-up of the polycrystalline structure (LD: loading direction).

Fig. 7. (a) ECCI image taken in the centre of the ½9 10 0-oriented grain after 100 cycles; ð111Þ and ð111Þ slip plane traces are visible. (b) Simulation set-up of the polycrystalline structure (LD: loading direction).

structure and by the boundary conditions, and can be calculated using this computational method. Dislocation structures are analysed in the central grain where the mesh is refined. In the simulations a pure tension-compression cyclic deformation, with a period of 0.004 s and a displacement amplitude of 1.48 mm (Fig. 5 (b)), 1.43 mm (Fig. 6 (b)) and 1.99 mm (Fig. 7 (b)), is applied along the loading direction (LD). The deformation amplitude is chosen such that the local strain amplitude in the central grain is 0:95% for the ½3 29 1 and ½1 19 0-oriented grains, while it is 0:65% for the ½9 10 0-oriented grain, as measured by in situ DIC. The strain rate is 9.5 s1 and 6.5 s1 for the two different strain amplitudes. The strain rate in the simulation is higher than the experimental one in order to reduce the computational time. In a strain controlled test, the distance travelled by dislocations depends only on the strain increment if the dislocation density and stress distribution are uniform, as can be found by integrating over time Orowan's law (3). When dislocation non-uniformities are present, the dislocation configuration can be affected by the strain rate if the effective shear stress taeff , defined in (8), becomes comparable with the difference between the threshold stress tath in low and high dislocation density regions. This is due to the additional stress necessary to move dislocations in low dislocation density regions, compared with the threshold stress tath , during a higher strain rate simulation. This higher stress can move dislocations in high dislocation density regions at a strain amplitude smaller than in lower strain rate simulations. The typical effective resolved shear stress in our simulations is around 20 MPa, while the difference in threshold stress is around 100 MPa. Additionally, experimental investigations made at a deformation frequency around 1 Hz Laird et al. (1986) and around 0.01 Hz Pham et al. (2011) show a similar geometry of the labyrinth structures. It is worth to note that our strain rate is smaller than the one used in discrete dislocation dynamics simulations, which is typically around 104 s1 (Fivel, 2008). Average values of initial dislocation densities in the simulations are summarized in Table 5 and random non-uniformities of 10% are added. The value of these non-uniformities depends on the position and is the same for every dislocation density of the models in order to obey connectivity. The initial curved dislocation densities are chosen smaller than the densities with zero curvature. This is an approximation of the equilibrium state of

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

113

Table 5 Initial dislocation densities in the simulations using the edge-screw and dislocation junction model. Edge-screw model

r⊥þ ¼ 1,1013 m2 rNþ ¼ 1,1013 m2 r⊥þ;Nþ ¼ 1,1011 m2 r⊥;Nþ ¼ 1,1011 m2

r⊥ ¼ 1,1013 m2 rN ¼ 1,1013 m2 r⊥þ;N ¼ 1,1011 m2 r⊥;N ¼ 1,1011 m2

Dislocation junction model

r==þ ¼ 1,1013 m2 r£þ ¼ 1,1013 m2 r==þ;£þ ¼ 1,1011 m2 r==;£þ ¼ 1,1011 m2

r== ¼ 1,1013 m2 r£ ¼ 1,1013 m2 r==þ;£ ¼ 1,1011 m2 r==;£ ¼ 1,1011 m2

this model, for which dislocation creation and annihilation are balanced. The same initial conditions have been used to model dislocation structures and mechanical hardening during single slip deformation (Grilli et al., 2015). If the ratio between the initial curved densities and the initial non-curved densities is higher, the hardening would take place in a few cycles due to the higher dislocation multiplication rate, being inconsistent with our experimental observations. Because of the chosen value of the curved dislocation densities and the parameter L ¼ 0:1 mm, the initial condition can be thought as a distribution of dislocation loops with radius 0.1 mm and density 1,1011 m2, which are typical values in a pristine metal (Sandfeld and Zaiser, 2015). The 6 mm thickness of the representative volumes is chosen based on the typical dislocation travel distance. If the initial dislocation densities in Table 5 are used in Orowan's law (3), then the dislocation travel distance is around 1 mm and the surface effects due to dislocations exiting the geometry through the bottom surface of the representative volumes in Figs. 5 (b) and 6 (b) are minimized. The dislocation velocity, according to the imposed strain rate and Orowan's law (3), is lower than 103 mm/s, which corresponds to the limit imposed on the dislocation velocity. The dislocation structures represent a barrier for the motion of dislocations because, inside high dislocation density regions, the dislocation velocity decreases. Thus, the velocity limit does not affect the behaviour of dislocations, apart from those in elements where the dislocation density becomes temporary close to zero. 4.2.1. Dislocation patterning in the ½3 29 1-oriented grain The slip traces visible in Fig. 5 (a) belong to ð111Þ and ð111Þ planes. As revealed by Schmid factor analysis, the two most active slip systems are ½011ð111Þ and ½011ð111Þ. This last slip system is taken as reference to define the dislocation junction ! vector l lock. The simulated total dislocation density after 50 cycles in the central grain is shown in Fig. 8 (a) using the dislocation junction model and in Fig. 9 (a) using the edge-screw model. Dislocation walls separated by low dislocation density regions are observed using both models. This density has been found by averaging the total dislocation density at a depth of 100 nm and 300 nm, i.e. in the centre of the two elements that are closest to the upper surface, in order to be comparable to the analysis depth of ECCI experiments. Indeed, the penetration depth of 30 keV electrons, used in the experiments, is around 400 nm Drouin et al. (2007) and the dislocation contrast decays exponentially with depth (Zaefferer and Elhami, 2014).

Fig. 8. (a) Total dislocation density in the centre of the ½3 29 1-oriented grain after 50 cycles using the dislocation junction model. (b) Bitmap conversion used to calculate the 2D orientation distribution function and dislocation line orientation angle q.

114

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

Fig. 9. (a) Total dislocation density in the centre of the ½3 29 1-oriented grain after 50 cycles using the edge-screw dislocation model. (b) Bitmap conversion used to calculate the 2D orientation distribution function and dislocation line orientation angle q.

Fig. 10. (a) Total dislocation density in the centre of the ½3 29 1-oriented grain using the dislocation junction model at a depth of 300 nm after (a) 50 cycles and (b) 100 cycles.

After 50 cycles, the simulated dislocation pattern is stable and only the maximum value of the dislocation density increases. This is shown in Fig. 10, where the total dislocation densities at 50 and 100 cycles are compared. The maximum dislocation density within the dislocation structures rpeak for all the slip systems is shown in Fig. 11 (a) for the edge-screw model and in Fig. 11 (b) for the dislocation junction model. This quantity represents the contribution of the different slip systems to the formation of dislocation structures. The slip systems ½011ð111Þ and ½110ð111Þ have higher rpeak for the edge-screw model, while the slip systems ½110ð111Þ and ½011ð111Þ give the major contribution for the dislocation junction model. In Fig. 11 the definition of the dislocation line orientation angle q for edge, screw, == and £ dislocations is reported. A quantitative analysis of the orientation of walls is possible using the 2D orientation distribution function (Gasparyan and Ohanyan, 2015). A ridge filter is applied to the simulated dislocation density field and continuous dislocation walls are obtained. Then the images are converted into bitmaps, in which every point can be either white (representing a dislocation wall) or black (representing a channel), as shown in Fig. 8 (b) and Fig. 9 (b). The 2D orientation distribution function counts the number of line segments of length L, oriented at an angle of q to the horizontal axis, for which 90% of the length lies in the white area, representing dislocation structures. The 2D orientation distribution as a function of the angle q is shown in Fig. 12 (b) for the dislocation junction model (continuous line) and for the edge-screw model (dashed line). A fixed length of about L ¼ 1 mm is used. The choice of L in an interval between 1 mm and 1.4 mm does not affect the position of the maxima of the 2D orientation distribution function. The 2D orientation distribution function has the two highest peaks around 40 and 60 for the dislocation junction model, orientations close to £ segments of the two slip systems with highest rpeak ½110ð111Þ and ½011ð111Þ, and around 50 and 135 for the edge-screw model, orientations close to edge and screw segments of the slip system with the second highest rpeak ½110ð111Þ. The other peak around 117 corresponds to an orientation close to edge segments of the slip system with highest rpeak ½011ð111Þ. The line segments oriented at q ¼ 45o for the dislocation junction model and at q ¼ 117o for the edge-screw

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

115

Fig. 11. Maximum dislocation density rpeak in the centre of the ½3 29 1-oriented grain after 50 cycles on all the slip systems using (a) the edge-screw model and (b) the dislocation junction model.

Fig. 12. (a) ECCI image in the centre of the ½3 29 1-oriented grain processed using a Gaussian filter to remove slip lines, dislocation structures have a dark contrast. (b) 2D orientation distribution function for the edge-screw and the dislocation junction model.

model are shown in Fig. 12 (b). In the filtered ECCI image in Fig. 12 (a) dislocation walls are mainly parallel to the ½110 direction (q ¼ 45o ) and no walls parallel to ½110 are evident, in better agreement with the orientation distribution function of the dislocation junction model. 4.2.2. Dislocation patterning in the ½1 19 0-oriented grain The slip traces visible in Fig. 6 (a) belong to ð111Þ, ð111Þ and ð111Þ planes. By contrast, the two slip systems with the highest Schmid factor (0.429) are ½011ð111Þ and ½011ð111Þ. The two slip systems with the second highest Schmid factor (0.406) ½110ð111Þ and ½110ð111Þ are likely to be active. This last slip system is taken as reference to define the dislocation ! junction vector l lock. The simulated total dislocation density after 100 cycles in the central grain is shown in Fig. 13 (a) using the dislocation junction model. subsec:resultsPatterning The maximum dislocation density rpeak is shown in Fig. 13 (b) and it is highest for the slip systems ½110ð111Þ and ½110ð111Þ. Thus, CPFE simulations can predict the activation of these two collinear slip systems and show that the influence of neighbouring grains is changing the slip activity in the central one with respect to an equally oriented single crystal. When the two collinear slip systems are active, the dislocation junction model and the edge-screw model are identical because screw and == dislocation lines coincide. The 2D orientation distribution function for the ½1 19 0-oriented grain is shown in Fig. 14 (b), where a length of about L ¼ 1 mm is used. Three maxima around 0o , 63o , 135o are present. The q ¼ 135o orientation corresponds to == dislocations of the two active collinear slip systems and q ¼ 63o to £ dislocations on the ð111Þ slip plane. In the filtered ECCI image in Fig. 14 (a) several orientations of the dislocation

116

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

Fig. 13. (a) Total dislocation density in the centre of the ½1 19 0-oriented grain after 100 cycles using the dislocation junction model. (b) Maximum dislocation density rpeak in the centre of the ½1 19 0-oriented grain after 100 cycles on all the slip systems using the dislocation junction model.

Fig. 14. (a) ECCI image in the centre of the ½1 19 0-oriented grain processed using a Gaussian filter to remove slip lines, dislocation structures have a dark contrast, examples of oriented walls are enclosed by rectangles. (b) 2D orientation distribution function for the simulated dislocation density in Fig. 13 (a).

walls are visible: 0o coincides with the orientation of £ dislocations of the ½011ð111Þ slip system, 27o and 90o with £ and == dislocations on the ð111Þ slip plane. Some other walls with an orientation around 170o can also be identified. 4.2.3. Dislocation patterning in the ½9 10 0-oriented grain The most evident slip traces visible in Fig. 7 (a) belong to ð111Þ and ð111Þ planes. As revealed by Schmid factor analysis, the most active slip systems are ½011ð111Þ and ½011ð111Þ, followed by ½101ð111Þ and ½101ð111Þ. The simulated total dislocation density after 8 cycles in the central grain is shown in Fig. 15 (a) using the edge-screw model. A dislocation pattern is developing and the maximum dislocation density within the dislocation structures rpeak for all the slip systems is shown in Fig. 15 (b). Dislocation walls at q ¼ 55o and q ¼ 125o are expected from dislocations belonging to the most active slip systems. These two peaks are visible in Fig. 16 (b) and correspond to the projections on the observation plane of dislocation walls perpendicular to the Burgers vectors of ½011ð111Þ and ½011ð111Þ. However, in the filtered experimental image in Fig. 16 (a) the most evident and long dislocation walls are parallel to the loading direction (LD). 4.3. Mechanical properties In order to compare the mechanical properties predicted by the edge-screw model and dislocation junction model with the experiment, we have analysed the maximum value of syy as a function of the number of cycles for the ½3 29 1- and the ½1 19 0-oriented grain. This component is the stress along the tensile-compression direction in the experiment. It is averaged over the left boundary, perpendicular to the load axis, of the representative volumes in Fig. 5 (b) and Fig. 6 (b). The results are

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

117

Fig. 15. (a) Total dislocation density in the centre of the ½9 10 0-oriented grain after 8 cycles using the edge-screw model. (b) Maximum dislocation density rpeak in the centre of the ½9 10 0-oriented grain after 8 cycles on all the slip systems using the edge-screw model.

Fig. 16. (a) ECCI image in the centre of the ½9 10 0-oriented grain processed using a Gaussian filter to remove slip lines, dislocation structures have a dark contrast. (b) 2D orientation distribution function for the simulated dislocation density in Fig. 7 (a).

Fig. 17. Maximum stress as a function of the number of cycles (a) for the ½3 29 1-oriented grain and (b) for the ½1 19 0-oriented grain. The stress is shown for both the edge-screw dislocation model and the dislocation junction model. The stress component syy is averaged over the left boundary of the representative volume (Fig. 5 (b) and Fig. 6 (b)) and over the centre of the analysed grains.

118

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

shown in Fig. 17. Hardening takes place during the first 20 cycles, followed by a saturation of the maximum stress. The maximum stress increases slightly faster for the dislocation junction model and it has a slightly higher value for the edgescrew model. This is consistent with the higher dislocation density inside dislocation walls found using the edge-screw model (Figs. 9 (a) and 8 (a)). This reflects the higher threshold stress tath for the edge-screw model. We average syy also over the refined 200 nm elements in the centre of the analysed grains, corresponding to the regions where the dislocation patterning is studied (Fig. 8 (a) and Fig. 13 (a)). This stress is lower than the one averaged over the left boundary. The maximum stress value found experimentally (Fig. 2 in Nellessen et al. (2015) is around 350 MPa in the saturation regime and it is reached after approximately 20 cycles. The higher stress in the simulation is due in part to the higher strain rate used, as explained in section 4.2. At this higher strain rate, some dislocations can reach the maximum velocity vmax, as explained in section 2, and therefore, more strain is accommodated by the elastic part of the deformation gradient in (1). This increases the stress value and shows that the developed model is strain rate dependent. The stress component syy in the centre of the analysed grains is lower partly because of the grain orientation and partly because of the presence of dislocation structures. As explained in section 4.2, low dislocation density regions have lower threshold stress tath and plastic deformation can be accommodated at lower applied stress. This is not the case in neighbouring grains because the element size used is not fine enough to capture dislocation structures. Moreover, the maximum stress is determined by the dislocation density at saturation, which depends on the parameters for dislocation multiplication and annihilation. The uncertainty on these parameters for 316L stainless steel can partly explain the discrepancy between the simulated and measured stress. This is reasonable since the simulated dislocation density is typically higher than the one observed in the ECCI images shown in section 4.2. 5. Discussion The simulation of the ½3 29 1-oriented grain is used to distinguish between the predictions of the dislocation junction and the edge-screw model. When using the edge-screw model, walls parallel to the edge dislocation segments of the two slip systems with highest Schmid factor form (Figs. 11 (a) and 12 (b)), which are mostly constituted by edge dislocations. This is due to the high density of edge dislocations on the active slip systems after many cycles. Indeed, the annihilation distance for edge dislocations b d ⊥ is lower than for screw dislocations, as reported in Table 1. Edge dislocations move along their Burgers vector and interact, forming dislocation walls. Thus, the simulated channels tend to assume the shape of edge dislocation segments that moved away. By contrast, the dislocation junction model predicts that r== and r£ dislocation densities have the same order of magnitude after many cycles because their annihilation rate equations, like (17), involve both distances b d ⊥ and b d N . r£ dislocations of the slip systems ½110ð111Þ and ½011ð111Þ get away from forming low density regions along directions ! parallel to the dislocation junction vector, which in this case is l lock ¼ ½110. Thus, many dislocation walls form perpendicular ! o to l lock and this orientation is at around q ¼ 45 to the loading direction, in agreement with the experimental results in Fig. 12 (a) and by L'Espe'rance et al. (1986), who found wall orientations intermediate between f100g and f210g. The double pseudo-polygonisation (DPP) arrangement theory of Dickson et al. (1986). predicts energetically favourable wall orientations perpendicular to the directions that bisect the obtuse and acute angles between the edge dislocation lines on the two slip systems. According to Table 1 in Dickson et al. (1986) these directions are ½001 and ½210 for the ½3 29 1-oriented grain. This second type of walls is oriented at q ¼ 27o with the loading direction, similar to the value q ¼ 45o predicted using the dislocation junction model. Therefore, the experimental image in Fig. 12 (a) cannot be used to discriminate between the Dickson's theory and the dislocation junction model. However, the edge-screw model predicts q ¼ 117o and q ¼ 135o as main wall orientations, which are not observed experimentally. The formation mechanism of dislocation structures proposed by Li et al. (2011). does not predict walls perpendicular to the ½210 or ½110 direction. This theory predicts that the interaction of edge dislocations of the two mainly active slip systems creates walls perpendicular to the sum and difference of the two Burgers vectors, which in our case are ½001 and ½010. Walls perpendicular to the ½010 direction (q ¼ 90o orientation) are more common for the edge-screw model than for the dislocation junction model, as shown by the 2D orientation distribution function in Fig. 12 (b), but they are not observed in the experimental image in Fig. 12 (a). By contrast, our dislocation junction model predicts that dislocation walls are formed when £ dislocations move and interact forming immobile structures. This behaviour cannot be captured by models based solely on edge and screw dislocations. It is important to note that the new developed model does not assume that only Hirth junctions are present. As explained in section 3.1, the primary slip system is assumed as reference for the orientation of dislocation junctions. Therefore, all possible dislocation junctions forming through the interaction of dislocations of the primary and secondary slip systems are considered in the model, including Hirth junctions and the Lomer-Cottrell junctions. While the Hirth junctions may not be the only mechanism controlling the microstructure formation, it has been shown that it is sufficient as an explanation and therefore most probably plays a relevant role in this formation. The simulation of the ½1 19 0-oriented grain is used to examine the case of two active collinear slip systems on conjugate slip planes, in our case ½110ð111Þ and ½110ð111Þ, not considered by the DPP theory of Dickson et al. (1986). By applying the DPP theory criterion, the dislocation walls result perpendicular to ½001 and ½110. This second type of walls would be oriented at q ¼ 135o with the loading direction but they are not evident in the experimental image in Fig. 14 (a), being q ¼ 170o the closest visible orientation. The q ¼ 135o orientation corresponds to one of the three maxima of the 2D orientation distribution function in Fig. 14 (b) and these walls are constituted by = dislocations forming collinear junctions. However, the collinear junction involves an annihilation of dislocation segments along the intersection of the two slip planes ð111Þ and ð111Þ Madec

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

119

et al. (2003) and this can explain why this orientation is not visible in the experiment. Dislocation walls oriented at q ¼ 0o are predicted by our dislocation junction model as formed by dislocations on the ð111Þ slip plane, while the DPP theory criterion is not able to explain them. The q ¼ 90o orientation is explained by our model as dislocation junctions at the intersection between the ð111Þ and the ð111Þ slip planes. The predicted q ¼ 63o orientation corresponds to £ dislocations of the slip systems ½011ð111Þ and ½110ð111Þ. The observed orientation at around q ¼ 27o , as shown in Fig. 14 (a), coincides with the orientation of £ dislocations on the ð111Þ slip plane. The dominant q ¼ 0o orientation in the experimental image in Fig. 14 (a), compared with the 2D orientation distribution function in Fig. 14 (b), can be due to the lack of dislocation junctions between the slip planes ð111Þ and ð111Þ in the model, forming dislocation walls with q ¼ 0o orientation. Indeed the slip system ! ½110ð111Þ, and not ½110ð111Þ, is taken as reference to define the dislocation junction vector l lock. The activity of the ð111Þ slip system can also be underestimated in the simulation because the subsurface microstructure is not considered. The simulation of the ½9 10 0-oriented grain is used to examine the case of 〈011〉 oriented crystals. In this orientation veinchannel structures perpendicular to the Burgers vector of the most active slip system are expected (Li et al., 2011). However, the projections of the Burgers vector of the two most active slip systems ½011ð111Þ and ½011ð111Þ on the observation plane are at q ¼ 45o . Walls constituted of edge or screw dislocations belonging to these slip systems would appear at q ¼ 55o and q ¼ 125o . None of these orientations is observed in the experimental image in Fig. 16 (a). This indicates that the observed dislocation walls may be induced by dislocation junctions between the slip planes ð111Þ and ð111Þ. This shows another case where the comparison between experiments and the edge-screw model is poor. It should be noticed that our models do not fully reproduce the dislocation density continuity in the walls (Figs. 8(a), 9(a) and 13(a)). Our element size (200 nm) is comparable with the experimentally observed thickness of labyrinth walls, while the spacing between two neighbouring dislocation walls is around 1 mm in copper (Li et al., 2011). This can also explain why in some parts of the simulated dislocation structures, the continuity of the walls is missing. Reaction-diffusion approaches model the motion of dislocations using diffusion and flux terms in the constitutive equations. In the model by Pontes et al. (2006), the motion of dislocations on two slip systems takes place along the two orthogonal Burgers vectors respectively. Consequently low density regions form perpendicular to the Burgers vector directions and not, as observed experimentally, perpendicular to the dislocation junction line. The lack of dislocation junctions in the model, therefore, leads to incorrect predictions. Reaction-diffusion models are valuable to understand the characteristic length scale of dislocation structures, which is determined by the balance between dislocation interaction strength and dislocation mobility. According to equation (1) in Pontes et al. (2006), the distance between dislocation walls increases if the coefficient that determines the amount of immobile dislocations, called pinning rate in their model, or the homogeneous steady state dislocation density decreases. This applies also to our simulations: the small annihilation distance for edge dislocations b d ⊥ leads to a higher homogeneous steady state dislocation density in the edge-screw model than in the dislocation junction model. Therefore a larger length scale is predicted by the dislocation junction model. Additionally, a higher average dislocation density is predicted by the edge-screw model, as shown by the higher volume fraction of dislocation structures in Fig. 9 (a) than in Fig. 8 (a). The characteristic distance between dislocation walls is in the range 0:3÷0:6 mm for the edge-screw model in Fig. 9 (a) and 0:6÷0:8 mm for the dislocation junction model in Fig. 8 (a), which better agrees with the experimental image in Fig. 12 (a). This characteristic distance depends on the constitutive equations, while the initial dislocation distribution can affect the spatial position at which dislocation walls develop. If the strain amplitude is increased in our model, the volume fraction of dislocation structures increases Grilli et al. (2015) and the spacing between neighbouring dislocation walls decreases, in qualitative agreement with the similitude principle (Sauzay and Kubin, 2011). In the specific configuration of the polycrystal simulations presented in this article, the similitude principle can be partially hidden by the surface effect. This is because, in higher strain amplitude simulations, more dislocations exit the geometry and the average dislocation density is lower. Therefore, an exact inverse proportionality between the stress and the distance between dislocation walls in these polycrystal simulations is not expected. Finally, the higher strain rate used in the simulation and the uncertainty on the dislocation density at saturation lead to an overestimation of the maximum stress. However, the number of cycles to reach the saturation stress is correctly predicted. 6. Conclusions In this paper we present an approach to introduce the formation of dislocation junctions in a continuum simulation framework. A specific dislocation density is introduced, which takes into account dislocations oriented along the intersection of two conjugate slip planes, informing the model about the specific motion of dislocation segments perpendicular to the junction line. This new model predicts dislocation walls, arising from a random initial dislocation distribution, oriented parallel to the ½110 direction, which were not captured by existing reaction-diffusion models Pontes et al. (2006) and by low energy theories. The interpretation of electron channeling contrast imaging analyses of fatigued austenitic steel is improved, because our model identifies the interactions that cause the formation of this type of observed dislocation walls as the interactions between dislocations parallel and perpendicular to the dislocation junction line. The spacing between dislocation walls is shown to be correlated with the average dislocation density after many cycles and dependent on the mechanisms included in the simulations. The introduction of the dislocation junction mechanism is shown to be essential in continuum constitutive equations when describing fatigue on a sub-micrometer length scale.

120

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

Acknowledgements This work has been supported by the Swiss National Science Foundation, project no 138240, and by the Deutsche Forschungsgemeinschaft, project SA 2251/2-1. We thank Dr. Franz Roters and Prof. Philip Eisenlohr for discussion on the CPFE method; Prof. Thomas Hochrainer and Prof. E. Aifantis for discussion on their respective models.

References Aoyagi, Y., Kobayashi, R., Kaji, Y., Shizawa, K., 2013. Modeling and simulation on ultrafine-graining based on multiscale crystal plasticity considering dislocation patterning. Int. J. Plasticity 47, 13e28. https://doi.org/10.1016/j.ijplas.2012.12.007. Arsenlis, A., Parks, D.M., 2002. Modeling the evolution of crystallographic dislocation density in crystal plasticity. J. Mech. Phys. Solids 50 (9), 1979e2009. https://doi.org/10.1016/S0022-5096(01)00134-X. Bitzek, E., Brandl, C., Derlet, P.M., Van Swygenhoven, H., 2008. Dislocation cross-slip in nanocrystalline fcc metals. Phys. Rev. Lett. 100, 235501. https://doi. org/10.1103/PhysRevLett.100.235501. Bonneville, J., Escaig, B., Martin, J., 1988. A study of cross-slip activation parameters in pure copper. Acta Metall. 36 (8), 1989e2002. https://doi.org/10.1016/ 0001-6160(88)90301-X. Catalao, S., Feaugas, X., Pilvin, P., Cabrillat, M.-T., 2005. Dipole heights in cyclically deformed polycrystalline {AISI} 316l stainless steel. Mater. Sci. Eng. A 400401 (0), 349e352. https://doi.org/10.1016/j.msea.2005.03.094. pre s, C., 2004. Mode lisation physique des stades pre curseurs de lendommagement en fatigue dans lacier inoxydable auste nitique 316l. Institut national De polytechnique de Grenoble. Ph.D. thesis. http://www.theses.fr/2004INPG0133. pre s, C., Fivel, M., Tabourot, L., 2008. A dislocation-based model for low-amplitude fatigue behaviour of face-centred cubic single crystals. Scr. Mater. 58 De (12), 1086e1089. https://doi.org/10.1016/j.scriptamat.2008.02.027. Devincre, B., Hoc, T., Kubin, L.P., 2005. Collinear interactions of dislocations and slip systems. Mater. Sci. Eng. A 400401, 182e185. https://doi.org/10.1016/j. msea.2005.02.071. Devincre, B., Kubin, L., Hoc, T., 2006. Physical analyses of crystal plasticity by DD simulations. Scr. Mater. 54 (5), 741e746. https://doi.org/10.1016/j.scriptamat.2005.10.066. Dickson, J., Boutin, J., L'Espe0 rance, G., 1986. An explanation of labyrinth walls in fatigued f.c.c. Met. Acta Metall. 34 (8), 1505e1514. https://doi.org/10.1016/ 0001-6160(86)90095-7. Drouin, D., Couture, A.R., Joly, D., Tastet, X., Aimez, V., Gauvin, R., 2007. CASINO v2.42-a fast and easy-to-use modeling tool for scanning electron microscopy and microanalysis users. Wiley Period. 29, 92e101. https://doi.org/10.1002/sca.20000. Fivel, M.C., 2008. Discrete dislocation dynamics: an important recent break-through in the modelling of dislocation collective behaviour. Comptes Rendus Phys. 9 (34), 427e436. https://doi.org/10.1016/j.crhy.2007.11.005. Franciosi, P., Zaoui, A., 1982. Multislip in f.c.c. crystals a theoretical approach compared with experimental data. Acta Metall. 30 (8), 1627e1637. https://doi. org/10.1016/0001-6160(82)90184-5. Gasparyan, A.G., Ohanyan, V.K., 2015. Orientation-dependent distribution of the length of a random segment and covariogram. J. Contemp. Math. Analysis 50 (2), 90e97. https://doi.org/10.3103/S1068362315020053. Grilli, N., Janssens, K.G., Swygenhoven, H.V., 2015. Crystal plasticity finite element modelling of low cycle fatigue in fcc metals. J. Mech. Phys. Solids 84, 424e435. https://doi.org/10.1016/j.jmps.2015.08.007. Hochrainer, T., Zaiser, M., Gumbsch, P., 2007. A three-dimensional continuum theory of dislocation systems: kinematics and mean-field formulation. Philos. Mag. 87 (8e9), 1261e1282. https://doi.org/10.1080/14786430600930218. Hochrainer, T., Sandfeld, S., Zaiser, M., Gumbsch, P., 2014. Continuum dislocation dynamics: towards a physical theory of crystal plasticity. J. Mech. Phys. Solids 63 (0), 167e178. https://doi.org/10.1016/j.jmps.2013.09.012. Jakobsen, B., Poulsen, H.F., Lienert, U., Almer, J., Shastri, S.D., Sørensen, H.O., Gundlach, C., Pantleon, W., 2006. Formation and subdivision of deformation structures during plastic deformation. Science 312 (5775), 889e892. https://doi.org/10.1126/science.1124141. Jin, N., 1983. Dislocation structures in fatigued copper single crystals oriented for double-slip. Philos. Mag. A 48 (5), L33eL38. https://doi.org/10.1080/ 01418618308236534. Jin, N., Winter, A., 1984. Cyclic deformation of copper single crystals oriented for double slip. Acta Metall. 32 (7), 989e995. https://doi.org/10.1016/00016160(84)90001-4. Jin, N., Winter, A., 1984. Dislocation structures in cyclically deformed [001] copper crystals. Acta Metall. 32 (8), 1173e1176. https://doi.org/10.1016/00016160(84)90123-8. € chen, K., Bo €hlke, T., Fritzen, F., 2008. On estimates for the effective shear modulus of cubic crystal aggregates. PAMM 8 (1), 10551e10552. https://doi.org/ Jo 10.1002/pamm.200810551. Korsunsky, A.M., Hofmann, F., Abbey, B., Song, X., Belnoue, J.P., Mocuta, C., Dolbnya, I., 2012. Analysis of the internal structure and lattice (mis)orientation in individual grains of deformed CP nickel polycrystals by synchrotron x-ray micro-diffraction and microscopy. Int. J. Fatigue 42, 1e13. https://doi.org/10. 1016/j.ijfatigue.2011.03.003 fatigue Damage of Structural Materials {VIII}. L. P. Kubin, Dislocation patterning, Mater. Sci. Technol.doi:10.1002/9783527603978.mst0051. Laird, C., Charsley, P., Mughrabi, H., 1986. Low energy dislocation structures produced by cyclic deformation. Mater. Sci. Eng. 81 (0), 433e450. https://doi. org/10.1016/0025-5416(86)90281-8 proceedings of the International Conference on Low Energy Dislocation Structures. Leung, H., Ngan, A., 2016. Dislocation-density function dynamics an all-dislocation, full-dynamics approach for modeling intensive dislocation structures. J. Mech. Phys. Solids 91, 172e203. https://doi.org/10.1016/j.jmps.2016.03.008. Leung, H., Leung, P., Cheng, B., Ngan, A., 2015. A new dislocation-density-function dynamics scheme for computational crystal plasticity by explicit consideration of dislocation elastic interactions. Int. J. Plasticity 67, 1e25. https://doi.org/10.1016/j.ijplas.2014.09.009. Li, X.-W., Zhou, Y., Guo, W.-W., Zhang, G.-P., 2009. Characterization of dislocation structures in copper single crystals using electron channelling contrast technique in sem. Cryst. Res. Technol. 44 (3), 315e321. https://doi.org/10.1002/crat.200800346. Li, P., Zhang, Z., Li, S., Wang, Z., 2009. Cyclic deformation and fatigue damage behaviors of ½1 4 14 oriented Ag single crystal. Philos. Mag. 89 (32), 2903e2920. https://doi.org/10.1080/14786430903130847. Li, P., Li, S., Wang, Z., Zhang, Z., 2011. Fundamental factors on formation mechanism of dislocation arrangements in cyclically deformed fcc single crystals. Prog. Mater. Sci. 56 (3), 328e377. https://doi.org/10.1016/j.pmatsci.2010.12.001. L'Espe0 rance, G., Vogt, J., Dickson, J., 1986. The identification of labyrinth wall orientations in cyclically deformed aisi-sae 316 stainless steel. Mater. Sci. Eng. 79 (2), 141e147. https://doi.org/10.1016/0025-5416(86)90397-6. Ma, A., Roters, F., 2004. A constitutive model for fcc single crystals based on dislocation densities and its application to uniaxial compression of aluminium single crystals. Acta Mater. 52 (12), 3603e3612. https://doi.org/10.1016/j.actamat.2004.04.012. Maaß, R., Petegem, S.V., Ma, D., Zimmermann, J., Grolimund, D., Roters, F., Swygenhoven, H.V., Raabe, D., 2009. Smaller is stronger: the effect of strain hardening. Acta Mater. 57 (20), 5996e6005. https://doi.org/10.1016/j.actamat.2009.08.024. Madec, R., Devincre, B., Kubin, L., 2002. Simulation of dislocation patterns in multislip. Scr. Mater. 47 (10), 689e695. https://doi.org/10.1016/S1359-6462(02) 00185-9.

N. Grilli et al. / International Journal of Plasticity 100 (2018) 104e121

121

Madec, R., Devincre, B., Kubin, L., Hoc, T., Rodney, D., 2003. The role of collinear interaction in dislocation-induced hardening. Science 301 (5641), 1879e1882. https://doi.org/10.1126/science.1085477. Manonukul, A., Dunne, F.P.E., 2004. Highe and lowecycle fatigue crack initiation using polycrystal plasticity, Proceedings of the Royal Society of London A: Mathematical. Phys. Eng. Sci. 460 (2047), 1881e1903. https://doi.org/10.1098/rspa.2003.1258. Martínez, E., Marian, J., Arsenlis, A., Victoria, M., Perlado, J., 2008. Atomistically informed dislocation dynamics in fcc crystals. J. Mech. Phys. Solids 56 (3), 869e895. https://doi.org/10.1016/j.jmps.2007.06.014. Messerschmidt, U., Bartsch, M., 2003. Generation of dislocations during plastic deformation. Mater. Chem. Phys. 81 (23), 518e523. https://doi.org/10.1016/ S0254-0584(03)00064-6. Misra, A., Prasad, R., Chauhan, V.S., Kumar, R., 2010. Effect of peierls stress on the electromagnetic radiation during yielding of metals. Mech. Mater. 42 (5), 505e521. https://doi.org/10.1016/j.mechmat.2010.01.005. Nellessen, J., Sandlbes, S., Raabe, D., 2015. Effects of strain amplitude, cycle number and orientation on low cycle fatigue microstructures in austenitic stainless steel studied by electron channelling contrast imaging. Acta Mater. 87 (0), 86e99. https://doi.org/10.1016/j.actamat.2014.12.024. Pham, M., Holdsworth, S., 2012. Dynamic strain ageing of AISI 316L during cyclic loading at 300 c: Mechanism, evolution, and its effects. Mater. Sci. Eng. A 556, 122e133. https://doi.org/10.1016/j.msea.2012.06.067. Pham, M., Solenthaler, C., Janssens, K., Holdsworth, S., 2011. Dislocation structure evolution and its effects on cyclic deformation response of AISI 316L stainless steel. Mater. Sci. Eng. A 528 (78), 3261e3269. https://doi.org/10.1016/j.msea.2011.01.015. Pham, M., Holdsworth, S., Janssens, K., Mazza, E., 2013. Cyclic deformation response of AISI 316L at room temperature: mechanical behaviour, microstructural evolution, physically-based evolutionary constitutive modelling. Int. J. Plasticity 47 (0), 143e164. https://doi.org/10.1016/j.ijplas.2013.01.017. Pirondi, A., Bonora, N., Steglich, D., Brocks, W., Hellmann, D., 2006. Simulation of failure under cyclic plastic loading by damage models. Int. J. Plasticity 22 (11), 2146e2170. https://doi.org/10.1016/j.ijplas.2006.03.007 damage and Fracture: Modeling and Experiments. Pontes, J., Walgraef, D., Aifantis, E., 2006. On dislocation patterning: multiple slip effects in the rate equation approach. Int. J. Plasticity 22 (8), 1486e1505. https://doi.org/10.1016/j.ijplas.2005.07.011 special issue in honour of Dr. Kirk Valanis Valanis Issue. Reed, R.P., Horiuchi, T., 1983. Austenitic Steels at Low Temperatures. Plenum Press. Roters, F., Eisenlohr, P., Hantcherli, L., Tjahjanto, D., Bieler, T., Raabe, D., 2010. Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: theory, experiments, applications. Acta Mater. 58 (4), 1152e1211. https://doi.org/10.1016/j. actamat.2009.10.058. Roters, F., Eisenlohr, P., Kords, C., Tjahjanto, D., Diehl, M., Raabe, D., 2012. Damask: the Düsseldorf Advanced Material Simulation Kit for studying crystal plasticity using an FE based or a spectral numerical solver. Procedia IUTAM 3 (0), 3e10. https://doi.org/10.1016/j.piutam.2012.03.001. Sandfeld, S., Zaiser, M., 2015. Pattern formation in a minimal model of continuum dislocation plasticity. Model. Simul. Mater. Sci. Eng. 23 (6), 065005. https://doi.org/10.1088/0965-0393/23/6/065005. Sauzay, M., 2009. Modelling of the evolution of micro-grain misorientations during creep of tempered martensite ferritic steels. Mater. Sci. Eng. A 510511 (0), 74e80. https://doi.org/10.1016/j.msea.2008.04.121, 11th International Conference of Creep and Fracture of Engineering Materials and Structures, {CREEP} 2008. Sauzay, M., Kubin, L., 2011. Scaling laws for dislocation microstructures in monotonic and cyclic deformation of fcc metals. Prog. Mater. Sci. 56 (6), 725e784. https://doi.org/10.1016/j.pmatsci.2011.01.006 festschrift Vaclav Vitek. Turner, A.P.L., Vreeland, T., 1970. The effect of stress and temperature on the velocity of dislocations in pure iron monocrystals. Acta metall. 18 (11), 1225e1235. https://doi.org/10.1016/0001-6160(70)90113-6. Vegge, T., Rasmussen, T., Leffers, T., Pedersen, O.B., Jacobsen, K.W., 2000. Determination of the of rate cross slip of screw dislocations. Phys. Rev. Lett. 85, 3866e3869. https://doi.org/10.1103/PhysRevLett.85.3866. Xia, S., El-Azab, A., 2015. Computational modelling of mesoscale dislocation patterning and plastic deformation of single crystals. Model. Simul. Mater. Sci. Eng. 23 (5), 055009. https://doi.org/10.1088/0965-0393/23/5/055009. Zaefferer, S., Elhami, N.-N., 2014. Theory and application of electron channelling contrast imaging under controlled diffraction conditions. Acta Mater. 75 (0), 20e50. https://doi.org/10.1016/j.actamat.2014.04.018.