Fourth-order strain-gradient phase mixture model for nanocrystalline fcc materials

The proposed modeling approach for nanocrystalline materials is an extension of the local phase mixture model introduced by Kim et al (2000 Acta Mater. 48 493–504). Local models cannot account for any non-uniformities or strain patterns, i.e. such models describe the behavior correctly only as long as it is homogeneous. In order to capture heterogeneities, the phase mixture model is augmented with gradient terms of higher order, namely second and fourth order. Different deformation mechanisms are assumed to operate in grain interior and grain boundaries concurrently. The deformation mechanism in grain boundaries is associated with diffusional mass transport along the boundaries, while in the grain interior dislocation glide as well as diffusion controlled mechanisms are considered. In particular, the mechanical response of nanostructured polycrystals is investigated. The model is capable of correctly predicting the transition of flow stress from Hall–Petch behavior in conventional grain size range to an inverse Hall–Petch relation in the nanocrystalline grain size range. The consideration of second- and fourth-order strain gradients allows non-uniformities within the strain field to represent strain patterns in combination with a regularization effect. Details of the numerical implementation are provided.

Keywords: nanocrystalline material, gradient plasticity, higher order gradient terms, constitutive modeling (Some figures may appear in colour only in the online journal)

Introduction
Nanocrystalline materials are defined as polycrystalline materials with an average grain size smaller than 100 nm. Their mechanical properties are distinctly different from those of coarse grained materials [2][3][4][5]. Therefore, bulk nanocrystalline materials receive great attention in the materials research community [6][7][8][9][10][11][12][13][14]. Of particular interest is the breakdown of the known grain size scaling relation, i.e. Hall-Petch relation [15,16], which holds for coarse-grained polycrystalline materials with grain sizes above µ 1 m, for which the strengthening effect is proportional to / d 1 , with d denoting the average grain size. A breakdown of the classical Hall-Petch relation for nanocrystalline materials was reported by several authors [3,11,13,[17][18][19][20][21][22], and the decrease of the yield strength with decreasing grain size is often referred to as an inverse Hall-Petch relation. This behavior is represented by a descending branch of the stress-/ d 1 -diagram. There have been several attempts and models to explain the origin of the inverse Hall-Petch relation, including grain/interface boundary sliding [13,23], diffusion over triple lines [24], absorption of dislocations by grain boundaries [25], grain-size dependent Frank-Reed dislocation sources [26], and composite models. Aifantis & Konstantinidis [27] provided a theoretical explanation of the inverse Hall-Petch behavior based on a gradient theory with an energy-related interface term. The same authors attributed the inverse Hall-Petch behavior to the presence of high interfacial area, and the existence of nanopores [28]. Trelewicz and Schuh [20] connected the breakdown to a change to inhomogeneous, glass-like flow behavior for very fine grain sizes. Argon and Yip [29] suggested that there are two competing mechanisms, grain boundary shear and dislocation glide, and associated the transition point to the inverse Hall-Petch behavior with the condition that both mechanisms contribute equally to the strain rate. In their model, the grain boundary shear mechanism is predominant below the transition point. Phase mixture models [1,24,30,31], treating the grain interior as one phase and the grain boundary as a separate amorphous second phase, successfully modeled and explained the transition from conventional to inverse Hall-Petch relation. As the main deformation mechanisms in nanocrystals, grain boundary diffusion and sliding were considered (e.g. [3,32,33]) which are closely related [34,35]. As stated by Kim & Estrin [35], it is essential to account for dislocation glide as well as diffusion mechanisms in models of the mechanical behavior of nanocrystalline materials. Upon further reduction of the grain size, the diffusion-controlled mechanisms are mainly responsible for plastic deformation. Bouaziz et al [36] pointed out that this effect led to a decrease in strain hardening and an increase in strain-rate sensitivity of the flow stress.
As nanocrystalline materials are highly sensitive to grain size and strain rate changes, inhomogeneities inevitably result in incompatibilities within the structure leading to gradient effects. Aifantis [37] showed that gradient enhancement (of second-order) for nanocrystalline materials captured the experimentally observed phenomena such as the size dependence of the elastic moduli and hardness. His results demonstrated a successful application of gradient terms in the context of elastic deformations, plastic flow, and diffusion.
Higher order gradients can be motivated by two considerations: (i) introduction of nonuniformities in the strain field to represent microstructural features more closely and (ii) smoothing of non-uniformities within the strain field to regularize the solution, see Askes et al [38]. Aspect (i) is represented by the second-order gradient model with a positive coefficient in the second gradient term whereas aspect (ii) is intended for incorporating the fourth-order gradient 5 . Fourth-order strain-gradient models have rarely been studied so far: Askes et al [38] proposed a model in the context of elasticity based on a discrete microstructure consisting of a chain of particles and springs. Estrin & Mühlhaus [48] introduced a fourth-order straingradient constitutive model for gradient plasticity. Both approaches result in an equivalent second-order gradient formulation but possess different features with respect to the fourthorder gradient. The solution of the displacement field with the (positive) fourth-order gradient elasticity model by Askes et al [38] exhibits harmonic terms multiplied with exponential terms whereas the solution of the model 6 by Estrin and Mühlhaus [48] only shows a superposition of harmonic terms.
In the following, we propose a second-and fourth-order strain gradient phase mixture model for nanocrystals. The model is implemented into an in-house finite element code to examine the capabilities of the chosen formulation. To the authors' knowledge, this is the first fourth-order strain-gradient phase mixture model for nanocrystalline materials presented in the literature.

Fourth-order strain-gradient phase mixture model for fine-grained materials
In the following, we extend the local phase mixture model by Kim et al [1] to include spatial non-uniformities and strain patterning. The point of departure is the stress-strain relation in a one-dimensional setting. Here, ε represents the local plastic strain and ( ) σ ε the local strain rate dependent stress. ∆ = ∇ div is the Laplacian operator; r and s are the coefficients of the second-and fourth-order derivative terms, respectively. In the following, a second-and fourthorder strain-gradient phase mixture model for fine-grained materials is presented and studied for various scenarios and boundary conditions. The polycrystalline material, an example of the microstructure is shown in figure 1, is subdivided into grain interiors and grain boundaries (see figure 2(a), assuming that the strain (and strain rate) is the same in both phases. The grain interiors are treated as one single phase. Regardless of the grain size d, the grain boundary width w = 1 nm is assumed to be constant [49,50]. Approximating the nanocrystalline grain as being cube-shaped [1], the volume fraction of the grain boundary phase is determined by 3 . Tschopp et al [51] investigated the volume fraction of grain boundaries and triple junctions for space-filling tetrakaidecahedral grains leading to a very similar distribution over the grain size as for cubic grains (see also [1]). The presented volume fraction is a simplification as only grain boundaries and grain interior are considered. 5 Alternatively, a negative sign of the second-order gradient term leads to a stabilized numerical behavior as well, which is usually the case in strain-gradient (crystal) plasticity models [37,[39][40][41][42][43][44] or non-local damage models [45,46]. In such models a fourth-order gradient (with a negative sign) would result in a destabilizing effect and therefore is omitted. A positive coefficient of the second-order gradient term leads to destabilizing behavior which motivates the incorporation of the fourth-order gradient (with a positive sign) to restore the well-posedness of the problem [47]. 6 For the case of gradient elasticity. Due to assumed parallel arrangement of grain boundary and the grain interior, the first stress contribution is additively determined by the rule of mixtures as ( ) where σ GB and σ GI denote the stress in the grain boundary and grain interior phase, respectively. The phase mixture model is augmented by gradient contributions as schematically illustrated in figure 2(b). Whereas the local phase mixture model by Kim et al [1] only allows for homogeneous solutions, the non-local gradient extensions capture non-uniformities.
Two gradient models are introduced, which are motivated by gradient models from the literature: a physical based fourth-order gradient phase mixture model (see e.g. [38,48]) as well as a postulated gradient enhancement of the model with a negative coefficient in the second-order gradient term (see e.g. [39][40][41][42]). Both gradient enhanced phase mixture models are able to capture the inverse Hall-Petch behavior, as will be shown later. Additionally, these extensions capture size-dependent strengthening and hardening. In the fourth-order gradient phase mixture model non-uniformities, e.g. in form of heterogeneities or boundary conditions, are introduced to represent the microstructure. The total stress is determined additively as  [ where h g is the gradient hardening modulus and the prefactor 2/5 is chosen following Askes et al [38]. The stress in the second-order strain-gradient phase mixture model reads: A drawback of the latter model is that it can only map very simple strain profiles. However, numerically, this model is less complex than the fourth-order gradient phase mixture model. The deformation mechanism in the grain boundary is assumed to be associated with diffusional mass transport along grain boundaries [1,35]. Thus, plastic flow in the grain boundary phase is determined by where Ω represents the atomic volume, k the Boltzmann constant, T the absolute temperature and A a material coefficient. Based on molecular-dynamics simulations, Yamakov et al [34] noted a dependence of A on the grain size, with values of up to A = 55.5 reported for the small grain size limit for spherical grain shape. A = 50 is assumed in this work. The coefficient of self-diffusion within the grain boundary is determined by Arrhenius' law: where Q GB denotes the activation energy comprising the energy for self-diffusion in the unstressed boundary and vacancy migration energy; D GB 0 is the preexponential factor in the grain boundary diffusivity [1]. The stress in the grain boundary phase controlled by diffusion cannot exceed a certain maximum value and thus considered to be capped at that level. Following Kim et al [52], this limiting stress equals the yield strength of the amorphous state σ am : with sgn • • • = | | ( ) / being the sign function. By this version of equation (6), the local model formulation by Kim et al [1] is extended to capture cyclic loading as well. Diffusion as well as grain boundary sliding are often considered as the main deformation mechanisms in nanocrystalline materials. As stated in [35], the grain boundary sliding mechanism requires grain boundary diffusion as an accommodation mechanism, so that both mechanisms are closely connected and the rate-controlling deformation mode is still grain boundary diffusion.
In the grain interior, dislocation glide mechanism (ε disl ) as well as diffusion mechanisms are considered to operate concurrently (following [1]), as the model is intended to cover several length scales, from coarse grained to nanometer grained materials. The relevant diffusion mechanisms are lattice diffusion, known as Nabarro-Herring creep [53], contributing the strain rate ε −N H , and the grain boundary diffusion. The latter is referred to as Coble creep [54], contributing the strain rate ε Coble . Consequently, the total plastic strain rate for the grain interior phase is given aṡ˙˙˙.
The dislocation glide contribution is determined by where ε 0 is a reference strain rate, σ 0 is the initial yield stress associated with the initial dislocation density ρ 0 , and 1/m describes the strain rate sensitivity of the flow stress. ( ) H • denotes the Heaviside step function, accounting for the limiting grain size d c for the dislocation-related viscoplastic model to be applicable. The evolution equation of the dislocation density ρ is given by [55] The quantity accounts for dislocation storage at grain boundaries, where M represents the Taylor factor, b the magnitude of the Burgers vector, α a material constant, and µ GI the shear modulus of the grain interior. The first and second terms in equation (9) account for strain hardening with C 1 being a constant. The last term describes dynamic recovery where C 20 denotes a material parameter, ε * a reference strain rate and n the (temperature-dependent) exponent associated with dynamic recovery. The strain rate due to Coble creep reads [1] k The strain rate due to diffusion within the bulk (Nabarro-Herring creep) is written as The set of equations (6)-(12) provide a full constitutive description of the local phase mixture model which is used with the gradient extended stress relations (3) or (4) to complete the higher-order strain-gradient phase mixture model.

Numerical treatment of gradient enhanced phase mixture model
To solve for the second-order gradient term in the finite element context, we apply a dual mixed approach in spirit of Svedberg & Runesson [56], where we introduce a global variable g 2 , denoting the strain gradient as .
For the implementation of the fourth-order gradient model, we introduce an additional global variable g 4 , representing the fourth gradient of the displacement following the implementation of the second-order gradient model. The weak forms of the boundary value problem involving the fourth-order gradient read which satisfy the boundary conditions for all test functions * u , * g 2 and * g 4 . L 0 denotes the length of a one-dimensional specimen. Equation (15) 1 represents the linear momentum balance with the stress σ given by equations (3) and (4), respectively. Equations (15) 2 and (15) 3 are the weak forms of the field equations for the strain gradients g 2 and g 4 , respectively. This approach allows to prescribe boundary conditions up to the third gradient of the strain ( ε = ∂ g x 4 3 ) and it is sufficient to work with linear interpolation functions. In turn, this results in a constant strain ε, second-order strain-gradient ε ∂ x 2 and fourth-order strain-gradient ε ∂ x 4 in each element. A proof of concept of the implementation procedure for gradient elasticity is presented in the appendix. It is useful to work with appropriately scaled forms of equation (15) within the finite element implementation, e.g. normalizing the implementation with respect to an average element size, as otherwise for small element sizes the stiffness matrix is ill-defined. This is especially true for the applied monolithic solution approach in this work. The gradient augmented phase mixture model involves the numerical solution for the total plastic strain rate, equation (7), which is challenging as it includes a power-law function via equation (8). As the model is based on the strain rate (as opposed to the strain itself), even a small time step can result in drastic changes of the strain rate, especially for large values of m. Consequently, the value from the last time step may not represent an adequate starting point for the local Newton iteration even if the time step is small. Following the method presented by Wulfinghoff and Böhlke [57], the power law form is approximated as linear for ε β ε >ḋ isl where β ε denotes a chosen strain rate value. Therefore, the strain rate is determined by a piecewise definition as GI,˙ denotes the stress at the strain rate β ε 7 . This approach is illustrated schematically in figure 3.

Numerical examples of higher-order strain-gradient phase mixture models
The numerical examples are based on a one-dimensional copper rod (length L 0 , cross-sectional area A 0 ) loaded in tension at different strain rates. The material parameters of copperthe material of choice in this study-are summarized in table 1. The rod is discretized into equally spaced elements. Linear shape functions are used for all global field equations. The left end of the bar ( | = = u 0 x 0 ) is kept fixed and a constant strain rate is applied to the right end 7 Considering that ˙ḋ isl ⩽ ε ε holds, and that all contributions of the grain interior strain rates are in the same direction, it is convenient to work with 1 β > where β is close to one.
). In the following, we study two sets of boundary conditions for the higher-order terms. We refer to them as either 'micro-hard' or 'micro-free', in the spirit of boundary conditions in strain-gradient plasticity theories [58]. Micro-hard boundary conditions read where the latter two conditions are relevant only for the fourth-order gradient phase mixture model. In contrast, micro-free boundary conditions are specified as

Second-order strain-gradient phase mixture model
Local models cannot account for any non-uniformities or strain patterns, i.e. such models describe the behavior correctly only as long as the overall solution is homogeneous. As soon as heterogeneities or incompatibilities are present, this results in gradient effects which can be captured with a non-local extension of such models [48]. Additionally, sample size dependence has been reported for nanocrystalline materials [61,62] which requires the incorporation of a length scale into the model. Stress-strain curves for the phase mixture model augmented with a second-order gradient term with a negative coefficient according to equation (4) are displayed in figure 4. The results  (16) 1 , by a piecewise defined function to improve the convergence of the Newton scheme due to improved starting values for ε disl . Due to the strong coupling of ε disl and ρ˙, more than one 'linear' step may be necessary to be within the convergence radius of the Newton scheme.  ), Coble creep is the dominant deformation mechanism, whereas otherwise dislocation glide is dominant. The gradient contribution slightly depends on the predominant deformation mechanism as this influences the resulting strain profiles.
reflect the size dependence of the stress which strongly depends on the internal length scale l. For l = 0, the local phase mixture model is recovered and the results qualitatively agree with the ones by Kim et al [1]. With increasing length scale l the stress level rises, in accordance with strain-gradient plasticity models in the literature [40,[63][64][65]. This effect is similar for all grain sizes as well as strain rates, the exact value of the gradient stress contribution depending on the local strain profile, see for example figure 5(a). Due to the micro-hard boundary conditions, which constrain the plastic strain, a heterogeneous deformation profile occurs over the length of the rod. The resulting boundary layer increases its width with increasing internal length scale l, which is consistent with the literature [40].
In the classical Hall-Petch region (here exemplarily seen for d = 100 nm and above), an increase of the grain size leads to widening of the boundary layer as well (see figure 5(a)). It should be noted that the predicted exponent in the Hall-Petch relation (d −n ) lies in the range of 0.6-0.7, i.e. is somewhat different from the canonical Hall-Petch value of 0.5. As known from previous studies on the local phase mixture model (e.g. [67]), for very fine-grained mat erials at low strain rates (here illustrated for ≈ d 50 nm and ε ≈ − − 10 s 5 1 ) an inverse Hall-Petch behavior is observed, which is captured by the higher-order strain-gradient model as well, as shown in figure 5(b). The numerical results show the 'strongest size', the transition grain size from classical Hall-Petch to an inverse behavior, in the grain size range d = 18-70 nm depending on the applied strain rate. The experimental results of [66] identify the strongest size as d = 22 nm. The inverse Hall-Petch effect leads to a diffuse boundary layer in the strain profile, where the sharpest strain profile with the smallest boundary layer is observed at the transition point. A further decrease of the grain size induces an increasing boundary layer thickness. An inverse Hall-Petch behavior is characteristic for nanocrytalline materials [3,19,66,68]. Two distinct types of behavior are seen. For coarse-grained material, the slope . An increase in strain rate leads to a shift to smaller grain sizes for the transition point between the two trends.
in the Hall-Petch diagram is positive while below a critical grain size the Hall-Petch relation breaks down and the slope becomes negative. The critical grain size is strongly strain rate dependent. As explained in Kim and Estrin [67], diffusion-controlled mechanisms are the main cause of the inverse Hall-Petch effect as, at a fixed strain rate, stress increases with increasing grain size according to equations (11) and (12). In contrast, for the dislocation mechanism the stress decreases with grain size as the rate of dislocation density accumulation is inversely proportional to the grain size as seen from equations (9) and (10). Strain hardening is accounted for by the evolving dislocation density ρ. As soon as the grain size approaches a (small) critical magnitude ( → d d c ) as well as for very low strain rates the dislocation contribution decreases and becomes negligible. Therefore, for these conditions the model behavior is completely controlled by diffusion mechanisms. In this limit, strain hardening vanishes and the stress level is constant and proportional to the strain rate. This implies an increasing strain rate sensitivity of the flow stress in nanocrystalline materials compared to their coarse-grained counterparts, a trend consistent with experimental observations for fcc materials (see, e.g. [69][70][71][72]).
In summary, the second-order strain-gradient phase mixture model captures size-dependent stress-strain behavior as well as the inverse Hall-Petch effect. However, the model is not capable of explicitly incorporating strain patterns. This drawback is overcome by the fourthorder strain-gradient phase mixture model considered below.

Fourth-order strain-gradient phase mixture model
In the following, numerical results for the fourth-order strain-gradient phase mixture model are discussed. The incorporation of a second-order gradient term with a positive coefficient provides the possibility to introduce non-uniformities within the structure, e.g. due to strain incompatibilities, to represent strain patterns accurately within the model. In turn, the fourthorder gradient regularizes the model leading to a stable numerical scheme.
As the influences of the gradient hardening modulus h g and length ratio l/L 0 are more complex in the fourth-order strain-gradient phase mixture model due to the interplay between both gradient contributions, these are investigated in more detail here. Figure 6 shows the numerical results for the variation of the gradient hardening modulus h g for micro-hard boundary conditions. As an example, a grain size of d = 10 nm and a strain rate of ε = − − 10 s  7(a)), which can be attributed to spontaneously evolved strain patterns. Strengthening effects in nanocrystalline materials due to sample size have experimentally been found in micropillar tests and a variation of sample size for fixed grain sizes (e.g. [61,62] . The strain profile (see figure 7(b)) shows two maxima near the ends of the rod. With increasing length scale, these merge at the rod's center to one single maximum, a similar distribution as for the second-order gradient model. Further increase of l/L 0 does not change the strain profile. Hence, this microstructure represents a stable state.  The contributions of different stresses involved and the influence of the internal length scale l are illustrated in figures 8 and 9. In this example, a grain size d = 100 nm and a strain rate ε = − − 10 s 5 1 are assumed to be fixed, where the dislocation glide accounts for approximately 75% of the external strain ε = 0.1, see figure 9(a). For ε ≈ 0, the deformation is nearly completely dislocation driven, see figure 9(b). Figure 9(a) illustrates that the dislocation contrib ution accounts mainly for inhomogeneous strain. Similar observations were made for different strain rates, although the ratio of dislocation and diffusion controlled contributions varies. The observation that the dislocation part mainly accounts for the inhomogeneous strain is true even if diffusional creep is the dominant mechanism, see figures 9(a) and 10(f). With growing strain, the contribution of the diffusion controlled part increases-which becomes slightly less prominent with increasing l/L 0 , figure 9(b). For dislocation dominated deformation the gradient terms do not lead to a pronounced strengthening effect, see figure 8(a), in contrast to the diffusion controlled case ( figure 7(a)). The strain profile is similar to that in the diffusion controlled case, however, the strain maxima are more pronounced. As micro-hard boundary conditions are applied, the strain and therefore the strain rate are zero at the boundary. This implies that the stress in the crystalline and the grain boundary phases is zero at the specimen boundary as well, figures 8(c) and (d). Both stress contributions follow the strain profile, the undulations in the profile being more pronounced for σ GB as the boundary phase is only deformed by creep which is directly proportional to the total local strain rate. The higher-order boundary conditions enforcing ε ∂ | = ∂ 0 x B 2 so that the stress contribution σ ∇ 2 is zero at the boundary as well. To satisfy the momentum balance, σ ∇ 4 has to be the counterpart of the stress at the boundary, see figure 8(f). The stress due to the fourth-order strain gradient acts against the one of the second-order strain gradient, as expected. The high stress level for / ≈ l L 1 0 is related to the stress due to the fourth-order strain gradient alone, see figure 8(f). For d = 10 nm a lower stress level is obtained, indicating an inverse Hall-Petch behavior for this region. At this grain size, the dominant deformation mechanism is diffusion creep (figure 10(f)), where the strain rate associated with dislocation glide is three magnitudes smaller. The strain rate due to dislocation glide exhibits a distinct profile (figure 10(f)) which is visible for the overall strain profile (figure 10(c)) as well. This implies that the strain rate due to dislocation glide maps its strain profile even if the contribution is small (d = 10 nm), see figure 10(f). As evidenced by equation (6), the stress in the boundary phase σ GB linearly increases with the applied strain rate ε until σ am is reached. As the stress is proportional to d 2 , the critical strain rate at which σ am is reached, increases as d decreases. For = d 10 nm, the stress in the grain boundary is negligibly small, whereas for the other investigated grain sizes, at least partially over the rod, the amorphous strength limit of boundary phase σ am is reached. As soon as dislocation glide becomes the dominant mechanism (figure 10(f)) the stress in the grain interior σ GI increases with decreasing grain size, see figure 10(e), which is mainly related to the increase of the forest dislocation density ρ, figure 10(b). This represents the classical Hall-Petch behavior. The increase of the dislocation density ρ with increasing grain size is less pronounced for the fourth-order gradient model compared to results of the second-order gradient model. This implies that the gradient contribution to the stress is higher than for the case if the second-order gradient is present alone. The spatial distributions of the first-order strain gradient ε = ∂ g x 2 ( figure 10(d)) and third-order strain gradient ε = ∂ g show an opposite slope character, responsible for the stress contributions of the two gradient terms. . The strain-rate ε −N H is negligible. (a) Individual strain rate contribution ε disl (black) and ε Coble (red) over the rod at ε = 0.1. (b) Averaged (along the rod) strain rates ε disl (black) and ˙C oble εˆ (red) versus applied strain ε. Dislocation glide is the dominant deformation mechanism, which is mainly responsible for the strain distribution. With increasing strain, the contribution of diffusion controlled deformation mechanism increases up to 25% of the total strain rate. This effect becomes less pronounced with increasing length scale l, indicating that this microstructure favors dislocation glide.
As a final example, the situation of a heterogeneity within the rod is investigated in combination with micro-free boundary conditions. As a heterogeneity, a 5% smaller grain size is assumed for a domain length L 0 /10 in the middle region of the rod, which is interpreted as a result of a mechanical or heat treatment of the material. This heterogeneous domain within the rod leads to a strain jump which is smoothed by gradient terms, see figure 11(b). The strain jump is either positive or negative, depending on whether the material in the heterogeneous domain is softer or stiffer than the surrounding material. Which situation applies depends on the previously discussed change of the Hall-Petch relation. The global stress level increases with increasing grain size in case of an inverse Hall-Petch relation whereas the stress decreases otherwise. Therefore, the stress in the grain interior σ GI (figure 11(a)) shows a decrease or an increase in this region, which results in slightly less pronounced hardening at the macrostructural level. It is worth mentioning that the strain jump is less pronounced for the fourth-order strain-gradient model compared to the same simulation with the second-order strain-gradient model, showing a stronger regularization due to the gradient terms in case of the fourth-order straingradient phase mixture model. This indicates that strong heterogeneities are required to form a stable pronounced pattern.
This analysis of the fourth-order strain-gradient phase mixture model shows the possibility of introducing non-uniformities in the strain field to accurately map the occurring pattern in combination with a regularization effect. The fourth-order gradient leads to a regularization which results in a stable numerical scheme allowing further analysis of the model. As the material behavior predicted by the phase mixture model is highly sensitive to grain size and strain rate, heterogeneities result in strain incompatibilities, which require the application of a gradient theory as presented in this section.

Choice of length scale based on the grain size
The length scale l presents an intrinsic/characteristic material parameter which significantly influences the model behavior as seen in the previous examples. Experimental observations show deviations of the classical mechanical behavior as soon as specimen dimensions become comparable to the intrinsic length scale (e.g. [73][74][75]). However, one long standing debate in the context of gradient models is the identification of the internal length scale parameter l. Often, the length scale l is chosen to represent the material's microstructure. Kouznetsova et al [76] suggested that the internal length scale is given by / = l h 12 2 2 (h: RVE size) for representative volume elements (RVE) within a second-order computational homogenization procedure based on the FE 2 -method. Obviously, the choice of the RVE size plays a crucial role in this approach. Nix & Gao [77] correlated the length scale with microstructural features where l is defined as l = L 2 /b with L the spacing between dislocation obstacles and b the Burgers vector magnitude. Aifantis [78] related the length scale to the average dislocation cell size based on torsion and bending experiments. Voyiadjis and Abu Al-Rub [79] suggested a variable length scale rather than a fixed one. The authors derived a length scale function depending on the deformation as well as the resulting microstructural features. In the same line, Lapovok et al [80] proposed a physical definition of the length scale on the basis of the current dislocation cell size which is modified according to the deformation history. In the context of gradient models for metallic glasses, the length scale parameter is related to the fracture process zone size [81][82][83]. For an extensive overview of the identification and quantification of the length scale in gradient models the interested reader is referred to [84] where it is stated that a current trend is to relate the length scale parameters to the heterogeneity of the material.
Following Estrin and Mühlhaus [48] and Estrin et al [85], the relevant 'material elements' are identified with individual grains so that l has the meaning of the average grain size. We assume this to be valid for nanocrystalline polycrystals, as strain incompatibilities between grains are significant for the mechanical response. For this purpose three differently sized specimens ( { } µ = L 0.1, 1, 10 m 0 ) consisting of grains in the nanosized range, 10 nm to 100 nm, are investigated. In figure 12 the corresponding results are summarized for three different strain rates and micro-hard boundary conditions. The results indicate that the stress increases with decreasing specimen size L 0 as the grain size is large compared to the specimen size. In particular, this effect is visible for d = 100 nm and L 0 = 100 nm, corresponding to a monocrystalline specimen (L 0 /d = 1). For d = 10 nm, no significant stress increase is visible as the number of grains is ten (L 0 /d = 10) so that the maximum investigated ratio of the length scale and the specimen size is l/L 0 = 1/10. For certain material parameter combinations, even a slight softening effect is observed (e.g. for d = 10 nm, L 0 = 100 nm and . This strain profile is representative for the other strain rates as well. Strong size effects are visible for d = 100 nm and L 0 = 100 nm as the structure consists of a single grain. Influence of the gradient is reduced with decreasing grain size leading finally to the local stress-strain results. The strain profile shows an evolving pattern according to the grain size to specimen size ratio. In the limiting case (here exemplarily shown for d = {50, 100} nm and L 0 = 100 nm), a stable pattern with a single maximum is obtained. ε = − − 10 s 5 1 , see also figure 6(a). The exact amount of the hardening/softening contribution additionally depends on the gradient hardening modulus h g . The heterogeneous strain profile for ε = − − 10 s 1 1 , see figure 12(d), is representative for all strain rates. For a very small number of grains (i.e. five or less), a stable strain pattern with a distinct maximum is observed. As the number of grains increases, the maximum decreases until three local maxima are present (see also figure 7(b)). The maximum in the middle of the rod, for symmetry reasons, further decreases, whereas the other two maxima move towards the respective rod ends. For specimens that are large compared to the grain size, the local solution of a constant strain over the bar is recovered, as the enforced micro-hard boundary conditions affect only the grains at the specimen ends. The effect of micro-free boundary conditions with heterogeneities within the specimen can easily be seen from the results presented in figure 11, and no further discussion of it is included here for brevity.
The identification of the length scale with the grain size provides a direct association of finite element mesh with the microstructure. Heterogeneities between the grains occurring in the specimen can be directly incorporated. The results show that the effect of boundary conditions or heterogeneities diminish as soon as the specimen is 'significantly large' compared to the grain size. By contrast, if the number of grains is small, these factors play a crucial role in the deformation behavior. Additionally, the fourth-order strain-gradient phase mixture model is able to predict the transition from the classical to the inverse Hall-Petch correctly.

Summary
The approach taken to a constitutive description of nanostructured materials is based on the phase mixture model [1]. Since the phase mixture model is highly sensitive to grain size and strain rate, heterogeneities in the material result in incompatibilities which require the application of strain gradients.
Two different strain-gradient extensions of the phase mixture model are presented and the corresponding results of their numerical implementation have been discussed. The secondorder strain-gradient phase mixture model, with a negative coefficient in the second order strain-gradient contribution to the stress, shows that this model captures the sample sizedependent stress-strain behavior as well as the inverse Hall-Petch effect (grain-size-dependent behavior) correctly. The results indicate that a decrease of grain size leads to a reduction of the width of the boundary layer near a tensile specimen ends until the breakdown of the Hall-Petch relation. However, the second-order strain gradient model is not capable of incorporating strain patterns showing complicated strain distributions. This drawback has been overcome by the fourth-order strain-gradient phase mixture model, which includes second-and fourth-order strain-gradient contributions to the stress. Since the phase mixture model is highly sensitive to grain size and strain rate, heterogeneities result in strain incompatibilities. To handle these, the application of strain gradients is required. On the basis of the fourth-gradient approach, initial non-uniformities within the strain field have been shown to lead to characteristic strain patterns. The boundary conditions, heterogeneities, and size effects play a crucial role if the length scale approaches the dimensions of the specimen. A physically reasonable choice for the intrinsic length scale is its identification with the grain size. Mathematically, non-uniformities are induced by the positive second-order strain gradient, whereas the fourth-order strain gradient provides a regularization effect. The model is capable of correctly predicting the transition of flow stress from the Hall-Petch behavior in conventional grain size range to an inverse Hall-Petch relation for nanocrystalline materials. Additionally, an initial strengthening effect due to a spontaneously evolved strain pattern for smaller sample size is captured. The inhomogeneous strain profile is mainly accounted for by the dislocation glide contribution to plastic strain.
The developed fourth-order strain-gradient phase mixture model is able to interpret and explain various experimental phenomena, such as the inverse Hall-Petch relation, change in strain rate sensitivity, size-dependent strengthening and strain pattern formation reported for ultrafine grained metals [86], where patterning during severe plastic deformation was observed. It is argued that the included effects are of main importance if the sample size approaches the grain size in the nanocrystalline regime.
With the growing importance of gradient structures for industrial applications, particularly in products with surfaces nanostructured by surface treatment [87][88][89], it becomes necessary to have a computational modeling basis capable of handling such gradients. The present models provide such a tool, which has been demonstrated to be physically sound and computationally viable.  Figure A1 shows a comparison between the analytical solution (A.2) and the numerical solution obtained with the solution algorithm described in equation (15). Both solutions show an excellent agreement for the displacement u, the strain field ε and the strain gradient ε ∂ x . This holds true for further fields, e.g. ε = ∂ g x 4 3 , as well. A linear stability analysis of this model reveals that the incorporation of the fourth-order gradient does not unconditionally guarantee a uniqueness of the solution in case of constant stresses, a situation present in the phase mixture model for creep deformation alone. An appropriate choice for the internal length scale l leads to a stable and unique solution with the incorporation of the fourth-order gradient.
The analysis in the context of linear elasticity illustrates that the proposed finite element algorithm with linear shape functions is very well suitable to model the effect of higher (arbitrary) order strain gradients and provides a straightforward implementation scheme as applied to the phase mixture model.