The hardness and modulus of polycrystalline beryllium from nanoindentation

Nanoindentation was used to compare properties of four industrial beryllium grades with dif- ferent purity. An extremely high variation of hardness was observed in all samples which ob-scured differences between samples. Analysis of the nanoindentation data in combination with SEM/EBSD measurements demonstrated that the crystallographic orientation of the indented grain was the major source of the wide variation in hardness, which was 2.5 times higher when the indentation direction was close to the [0001] c-axis of beryllium compared to indentation along the [112¯0] or [11¯00] directions. The most noticeable difference between tested grades were observed for the “soft” orientations: hardness of less pure structural grades was 15–30% higher compared to the pure nuclear grade. Crystal plasticity finite-element (CPFEM) simulations in- dicated how the hardness anisotropy of beryllium arises from the anisotropy in the plastic deformation. Experiments and simulations also demonstrated that localised plastic deformation of the surface around the indent (pile-up or sink-in) was highly crystallographically dependent: during indentation into “soft” orientations, pile-up dominated increasing the contact area; while sink-in behaviour was dominant during indentation into “hard” orientation reducing the contact area. This implies that the hardness values calculated from indenter displacement and indenter profile using the standard Oliver-Pharr approach, without considering pile-up/sink-in effects, will be incorrect. Several contact area correction methods were applied and compared. In con- trast the indentation modulus was similar for all investigated grades and was not found to have any strong crystallographic dependence. Crystal plasticity finite element analysis indicates that this is due to the complex 3-dimensional nature of the elastic interaction between the indenter and the sample, and also since, for the chosen indentation depth, the elastic interaction volume is much larger than the materials’ grain size. surface orientation was independent of temperature. The lengths of the two diagonals of the Vickers impressions in (101¯0) and (112¯0) planes did change with temperature, which the authors linked to the alteration of pre-dominant slip systems. The results of Hill and Jones (1961) showed that work hardening of beryllium single crystals led to highly non-uniform increase of hardness: ∼30% for basal plane indentations but more than 3 times higher for indentations on prismatic planes. This paper adapted the approach for modern industrial beryllium with grain sizes too small for typical indentation tests used in macro- and micro-hardness measurements.


Introduction
Beryllium has a unique combination of mechanical and physical properties, including low density, low nuclear interaction crosssection, high strength and melting point combined with superior heat conductivity (Beryllium Science and Technology., 2012;Walsh, 2009) and it is successfully used in a wide variety of applications involving severe environmental conditions (Design Guides for Materion Beryllium Metal, 2017). Beryllium is of particular importance for nuclear application: it is one of very few candidate materials suitable for beam windows and target components in the next generation of neutrino sources (Hurh et al., 2013;Kuksenko et al., 2017), first wall material in the currently constructed International Thermonuclear Experimental Reactor (Mitteau et al., 2013) https://doi.org/10.1016/j.ijplas.2018.12.008 Received 26 September 2018; Received in revised form 12 December 2018; Accepted 18 December 2018 using the continuous stiffness measurement (CSM) technique (Oliver and Pharr, 1992) with a 42 Hz and 2 nm oscillation. The indenter tip was calibrated in fused silica to determine the contact area as a function of indenter displacement into the surface. The scanning mode of the G200 indenter was used to investigate the topography of the prints after indentation. The surface profiles were then extracted to correct the contact area due to pile-ups and sink-ins around the indent using open source software Gwyddion ("Gwyddion -Applications -http://gwyddion.net/," 2018).
To quantify nanoindentation hardness H and indentation modulus E, the method used as a first approach was the now-standard method proposed by Oliver and Pharr (1992) which forms the basis of the indentation testing standard (ISO 14577-1:2015. Metallic materials -Instrumented indentation test for hardness and materials parameters -Part 1: Test method, 2015). The following equations were used: where A is the projected contact area of the indentation at load P, v is Poisson's ratio of the test piece, 0.05 for beryllium (Beryllium Science and Technology, 2012), v i is Poisson's ratio of the indenter, 0.07 for diamond (ISO 14577-1:2015. Metallic materials -Instrumented indentation test for hardness and materials parameters -Part 1: Test method, 2015), E eff is the effective modulus of the contact between the indenter and the sample, E i is the modulus of the indenter, 1140 GPa for diamond (ISO 14577-1:2015. Metallic materials -Instrumented indentation test for hardness and materials parameters -Part 1: Test method, 2015), S is the contact stiffness, which is computed as the slope of the unloading curve continuously evaluated in the CSM mode, and is the tip geometry correction factor (1.034 for a Berkovich tip) (Fischer-Cripps, 2004).
In this study, the standard Oliver-Pharr algorithm for nanoindentation data treatment was modified in order to take into account the influence of the plastic deformation on the sample surface around indents during loading, by which the contact area is modified from that calculated via the standard method. Both the "Oliver-Pharr" and the corrected area values were used in equations (1) and (3) and the results are compared in section 3.2.5. In this paper H O-P , A O-P and E O-P refer to data obtained with the Oliver-Pharr algorithm, whereas H C , A C and E C -are values after the contact area correction.
At least 200 indents were made for each beryllium grade with a target maximum indentation depth of 400 nm. Surface position correction was made for each indentation before data analysis. For consistency, load and contact stiffness data averaged between 320 and 360 nm indentation depth were used for indentation hardness and modulus calculations, avoiding the depth range within which there is a strong indentation size effect. In order to investigate the crystallographic dependence of nanoindentation hardness, EBSD analyses of the indentation arrays were made, using a JEOL 840A scanning electron microscope equipped with the EDAX-TSL EBSD; OIM TSL software ("OIM Analysis™/EDAX," 2018) was used for analysis of the EBSD data. Surface normal-projected inverse pole figure orientation maps superimposed with image quality (a parameter quantifying the quality of the corresponding diffraction pattern (Wardle et al., 1994)) maps were used for indent localisation. Data from indents covering grain boundaries and in close proximity to them were analysed separately and were not considered in the crystallographic orientation analyses.

Crystal plasticity finite element method (CPFEM)
Finite element simulations using Abaqus 2016 ('Abaqus, 2016Documentation', 2016 were performed to investigate the influence of crystallographic orientation of the indented beryllium grain on the mechanical response. Simulations provided load as a function of indentation depth and details of the plastic deformation of the free surface around indents; deformation affects the contact area between the tip and the sample and hence the accuracy of the hardness and elastic modulus. A strain gradient crystal plasticity user material (UMAT) for Abaqus was developed based on the user element (UEL) written by Kiwanuka (Dunne et al., 2012, 2007). To allow the curl of the plastic deformation gradient to be calculated within a UMAT, the initial Gauss point coordinates and current value of the plastic deformation gradient F p at every integration point were stored in a common block; this was protected using the Abaqus utility subroutine MutexLock() to enable parallelisation. Full details of the curl computation are in Appendix A. A total of 193,800 elements (8 Gauss points/element) were used to represent a cube with dimension of L = 40 μm. A very fine biased mesh under the indent to be used: w = 50 nm, increasing up to w = 2 μm away from the indent. The indenter tip was modelled as a rigid object with a perfect Berkovich geometry. The small sliding, node to surface, Abaqus contact algorithm was used with the default frictionless hard contact property. The top (blue) surface was traction-free while the remaining 5 (red) surfaces were fixed as shown in Fig. 1 (c).
The total deformation F is decomposed multiplicatively into plastic (F p ) and elastic (F e ) deformation gradients: where the flow rule has the form The plastic velocity gradient, L p , is given by the crystallographic strain rate resulting from dislocation glide on the active slip systems with slip direction s and slip plane normal n .
The crystallographic slip rate on slip system α is given by (Dunne et al., 2012(Dunne et al., , 2007 where the prefactor A = 3 × 10 −3 /s for < a > slip and A = 1 × 10 −2 /s for < c+a > was treated here as a fitting parameter, although it can be related to the physically based constants , implying a mobile dislocation density ( m =3 × 10 −3 /μm 2 ), attempt frequency ( f =10 11 /s) Helmholtz free energy ( F = 3.46 × 10 −8 pJ) and Burgers vector length b = 0.228 nm for < a > slip and 0.424 nm for < c+a > . The activation volume was assumed to be = where B is a dimensionless fitting parameter (B = 0.05) and 0 an initial sessile dislocation density (assumed to be 1/μm 2 ) required to ensure a finite strain in the absence of geometrically-necessary dislocations (GND). Minimal fitting was performed: A and B were calibrated using one load displacement curve only; all other parameters were from literature. Hardening occurs through the activation volume parameter V (see Fig. 1(b)) which was assumed to be proportional to the total (edge + screw) GND spacing. A total of 12 slip systems was included in the model giving a possible 12 edge + 9 screw dislocation types. Only the total scalar density was used in the slip law which was obtained by summing the absolute values of all 21 densities, = = = GND 1

21
. The 21 densities were found by reformatting the dislocation density tensor (the transpose of the Nye-Kröner-Bilby tensor) into a set of linear equations where column α of M is the tensor product l b of the dislocation line direction and Burgers vector of slip system α, reshaped as a column vector, and the RHS is the curl of the plastic deformation gradient (Das et al., 2018).
also reshaped as a column vector (Arsenlis and Parks, 1999). Note that the reshaping of l b and G km into column vectors can be done row-wise or column-wise provided it is done consistently on both sides of equation (8). As M is not square (having 9 rows and 21 columns) it cannot be inverted; instead the pseudo right inverse was found which allows the L 2 minimum solution for the GND density on each slip system, , to be obtained for every possible slip system The resulting crystallographic slip rate is shown in Fig. 1 for different values of GND . It should be noted that twinning, which is known to be an important deformation mechanism in beryllium crystals under compression (Beryllium Science and Technology, 2012; Brown et al., 2012) was not incorporated in the current crystal plasticity model. Fig. 1. a) the < a > basal, < a > prismatic and secondary < c+a > slip systems. (b) The slip rate on the < a > basal systems for different geometrically necessary dislocation (GND) densities. As the GND density increases the slip rate decreases for a given stress due to GNDs acting as obstacles to dislocation slip. (c) The FE mesh.

Results and discussion
3.1. Microstructure Fig. 2 shows a comparison of the microstructures of the investigated beryllium grades. The white contrast objects on the micrograph are precipitates, typically enriched with oxygen, aluminium, iron, silicon, titanium and other impurities, as was demonstrated by our previous STEM/EDS studies of beryllium (Kuksenko et al., 2017). There is a noticeable difference between grades which is likely to have originated from the different thermomechanical treatment histories of the grades: in the hot rolled PF-60 grade precipitates are mainly located in small, often elongated agglomerates inside grains, whereas in the grades produced by VHP and HIP the precipitates are mainly at grain boundaries.
Grains mainly have irregular but roughly equiaxed shapes, containing some thin elongated twins. In this study the twins were not quantified, but from the visual inspection (see Fig. 2) it was noticed that the HIP S-200-FH grade has a larger number density of twinned grains, visible as thin elongated grains with inverted contrast inside equiaxed grains on Fig. 2(d). Fig. 3 shows EBSD maps of the investigated grades. Fig. 3(a) shows the area of the nanoindentation array in the PF-60 grade, where the normal-projected inverse pole figure orientation maps (in colour) is superimposed on the image quality maps (Wardle et al., 1994) (in grey). Dark areas indicate low EBSD image quality and correspond to grain boundaries and areas around nanoindentations. The SEM micrographs and EBSD maps demonstrate that each of the investigated grades has a considerable variation of grain size. PH-60 grade has noticeably larger grains than the other investigated grades. Average grain diameter measured from the EBSD maps is (25 ± 14) μm for the hot-rolled grade, (7.3 ± 4.4) μm in the S-65, (7.1 ± 4.1) in the S-200-F grade and (5.3 ± 3.3) μm for the HIP-ed S-200-FH grade. Note that since the step size for the EBSD mapping was from 0.4 to 0.6 μm, very thin twins were not considered during the grain-size calculations.
The EBSD inverse pole image in Fig. 3(a) shows that the PF-60 hot rolled beryllium grade has a crystallographic texture with a preferential orientation of the basal plane (0001) close to the plane of observation, as was also demonstrated in an earlier detailed microstructural study (Kuksenko et al., 2017). This texture results from the hot rolling used for forming the PF-60 sheet ("PF-60 Beryllium Foil for X-ray," 2018). The other investigated beryllium grades have no any noticeable texture.

Nanoindentation data
As anticipated, significant scatter in the hardness values was observed in all beryllium grades investigated ( Fig. 4(a)). The data exhibit an extremely broad distribution, with values ranging from 2.5 to 6.5 GPa. The shapes of the frequency distribution are different in the rolled and pressed grades: VHP and HIP grades have symmetric distributions with the maximum fraction at 3-5 GPa, whereas the hot rolled PF-60 distribution is skewed towards higher hardness values with a maximum fraction at 5-6 GPa. Average hardness values are 4.5 ± 0.8, 3.7 ± 0.6, 4.1 ± 0.5 and 4.5 ± 0.7 for the PF-60, S-65, S-200-F and S-200FH grades correspondently.   V. Kuksenko et al. International Journal of Plasticity 116 (2019) 62-80 4.5 GPa, while the S-65 grade has the lowest average hardness (3.7 GPa). Fig. 4(b) also compares the average harnesses from grains and grain boundaries and demonstrates a noticeable difference between these areas. The average hardness values from grain boundaries were 6-13% higher than at grain cores, that might originate from extra stresses needed to transfer dislocations across grain boundaries. Another possible source of the increased GB hardness is the impurity precipitates which are mainly located at grain boundaries in VHP and HIP-ed grades (see Fig. 2(b, c), and d). We should note that the hot rolled PF-60 grade the grains boundaries are free of precipitates ( Fig. 2(a)), but it is know, that in this grade grain boundaries have noticeable iron segregation (Kuksenko et al., 2017) which possibly can contribute to GB areas hardening. Our data (Fig. 4) demonstrated that both the average and distribution of hardness values vary between grades, however, data scattering and texture effects complicate the grades comparison. Therefore, a more detailed investigation of crystallographic dependence of indentation results was performed.

Load vs displacement data analysis
For crystallographic effects in indentation data analysis, loads as a function of indentation depth were grouped according to the angle between the indentation direction and the [0001] axis of the indented grain (which for the remainder of this paper is denoted ), as demonstrated in Fig. 5 for the case of the PF-60 beryllium grade. Indentations 400 nm deep into the basal plane (i.e. with small ) require about 17 mN, almost twice as high a load as for indentation into surface planes perpendicular to the [0001] axis (i.e.°9 0 ). The indentations made in three crystallographic directions were also simulated with CPFEM. The simulated load displacement curves for indentations into "hard" [0001] ( = 0) and "soft" [1010] and [1120] crystallographic directions = ( 90) (i.e. into surfaces normal to these directions) are shown in Fig. 5(b) as wide pink, green and dashed purple lines. The two free parameters in the slip law, the slip rate prefactor A and the hardening parameter B, relating GND density to activation volume, were calibrated to reproduce the load displacement curve in the [0001] orientation, with all other parameters taken from the literature. These parameters were then used for simulation of indentations in the [1010] and [1120] directions. The excellent agreement between experimental data and modelling for the two soft orientations implies that the CRSS values and slip law are reasonable and confirms that the origin of the experimentally observed scatter in nanoindentation results in these polycrystalline materials is anisotropy in the plastic slip from indentations into differently-orientated grains. Fig. 6(a)-(c) shows CPFEM contour plots of the total geometrically-necessary dislocation (GND) density, GND , on a log scale, at maximum load for indentations in the "hard" and the two "soft" orientations. The majority of GNDs are produced in a localized area just under the indent (yellow, orange and red zones on Fig. 6). In the hard orientation the region of high GND density is predicted to penetrate further into the sample than for the soft orientations. Easy < a > slip on basal and prismatic planes produces lobes of high plastic strain and hence high GND density in the 3 < a > directions and in the 3 directions made up of the combinations of two different < a > directions. This is shown in Fig. 6(c) where the soft orientation (indentation into a prismatic plane parallel to the free surface) produces a GND density pattern reflecting slip plane orientations when viewed down the c-axis. Fig. 6(d)-(f) shows the effective plastic strain at maximum load for simulations of indentations made in the 3 different crystallographic orientations. The plastic strain is highly localized near the indent, and therefore only a smaller region of material is shown than in Fig. 5(a)-(c). The maximum plastic strain value is greatest (> 40%) for the hard orientation directly under the indent. Field plots at different depths showed that the plastic zone depth and surface extent, defined by the volume with an effective plastic strain greater than 5%, was approximately 7 times larger than the tip penetration depth. This suggests that if the surface distance between the indent and the nearest sub-surface grain boundary is more than 8 times the indentation depth, then plastic response of the material will be that for a single crystal of the chosen crystallographic orientation. This approach was used in the current work to remove grain-boundary indentation data from single crystal data. Since it is not possible to determine grain structure in volumes under the surface cross-section, the approach naturally suffers from scattering generated by the possible presence of grains with different orientation in close proximity under the surface cross-section. Fig. 7 shows the predicted four in-plane components of the elastic stress tensor on the y = 0 plane directly under the indent for the hard [0001] and soft [1010] orientations at maximum load. Fig. 7(a) and (c) show that the soft orientation has a larger volume of material with a compressive xx , whereas the zz component of stress for the hard orientation penetrates very deep into the material compared to the soft orientation. The easy < a > slip in the soft orientation also reduces the elastic stress and hence reaction force, in good agreement with the measured load-displacement curves, as shown in Fig. 5 (b). The simulation results demonstrate that the field of significant elastic interaction between the indenter tip and beryllium samples is almost two orders of magnitude larger than the tip penetration depth and thus for the ∼400 nm deep indents used here it extends up to ∼40 μm from the indented surface. This is significantly larger than the average grain radius of the tested beryllium grades (from 2.7 to 12.5 μm).

Contact area
In order to calculate accurate hardness values, the contact area between the indenter and the sample in equation (1) should be correctly determined. In the standard method proposed by Oliver and Pharr (1992) the contact area used in equations (1) and (3) (Oliver and Pharr, 1992) where P max is the measured maximum load, S is the elastic unloading stiffness and is the Berkovich tip geometry coefficient. However it is known that elastic, plastic and strain-hardening properties of a material can affect its behaviour underneath the indenter tip and the material may either plastically sink-in or pile-up around the indenter (Bolshakov and Pharr, 1998;Fischer-Cripps, 2004;McElhaney et al., 1998). Both piling-up and sinking-in may significantly change the contact area between the indenter tip and surface (see examples from tungsten (Beck et al., 2017), titanium (Britton et al., 2010;Su et al., 2016;Zambaldi et al., 2012), magnesium (Haghshenas et al., 2018), nickel (Renner et al., 2016), Zircaloy (Kese et al., 2018) and ironbased alloys (Hardie et al., 2015)). If unaccounted for, the calculated indentation hardness and modulus would be too high if pile-up occurs (as the real contact area is greater than that calculated from the indenter profile and penetration depth) and too low when sink-in occurs (when the real contact area will be less than the calculated one).
For estimation of the out-of-plane plastic deformation of beryllium and hence to improve the analysis of the load-displacement data to give indentation hardness, direct imaging of the residual impressions, using indenter scanning under a low contact load (1 μN), was performed for 47 indents covering a range of indentation orientations for the 4 grades studied. Fig. 8 shows examples of experimentally measured topography maps around indents into 3 different grain orientations and CPFEM simulations for similar indentation orientations. For all indentations, plastically deformed surface zones were observed well outside the area of the tip penetration.
Estimations of the real contact areas were made using the approach proposed by Kese et al. (2005). The method requires experimental measurement of the indent profile through the central point of indentation and the centre of each edge, as shown in Fig. 8(a) and (b). Where pile-up is observed, Fig. 8(g), the contact between the indent and the pile-up is determined and is then projected onto the direction of the free surface normal, giving the pile-up height, h p . The contact area of each pile-up is then approximated as a semi-ellipse, with the major axis equal to the length of the side of the projected triangular area of the indent print, and the minor axis being measured on the indent profile image as the projected distance of the pile-up contact perimeter as demonstrated in Fig. 8(h). Originally, Kese et al. (2005) proposed the method for estimation of pile-up contribution only, however, in this work, the approach was modified to account for the observed sink-ins: the reduction in contact area was also approximated as a semi-ellipse, with a major axis calculated in the same way, whereas the minor axis was measured on the indent profile image as the projected distance between the sink-in contact perimeter as shown in Fig. 8(g). For each investigated indent the corrected projected area, A corr , was calculated as the contact area determined by the Oliver-Pharr analysis (Oliver and Pharr, 1992), A O-P plus (or minus) the area due to out-of-plane deformation near each edge for pile-ups (or sink-ins). Fig. 9(a) shows the magnitude of the out of plane plastic deformation for all the measured profiles, where 3 data-points across each edge of the investigated indent print were sorted from the smallest value (blue data points) to largest (red data points). The graphs show that, for all indentation directions, both pile-ups and sink-ins are simultaneously present in the vast majority of indentations. For indentation with < 15°(hard orientations), at least two edges had adjoined sink-ins whereas the third edge had either sink-in or small pile-up. The residual sink-in area was spread out up to ∼1.5 μm around the indent, as shown by the blue "halo" around the indenter print on Fig. 8(a) with the maximum depth up to ∼55 nm. The simulated residual out of plane displacement from the crystal plasticity model (Fig. 8(d)) qualitatively corresponds to the experimentally observed behaviour; however, the predicted deformed area was more symmetric than found in experiments. For indentations with ∼90°(soft orientations), pile-up dominated and exhibited two hills on opposite sides of the indentation pit. The pile-up maxima are aligned with the traces of the < a > directions onto the free surface. Depending on the orientation of the indenter tip two types of pattern were observed: i) If the c-axis of the grain was approximately perpendicular to one of the edges of the tip, this edge typically exhibited sink-in, and pile-up occurred at the other two edges (see Fig. 8(c)). ii) If one of the tip's edges was approximately parallel with the c-axis of the indented grain, then one large pile-up (up to 220 nm height) was observed near this edge, and the second pile-up was near the opposite corner of the indentation tip (see Fig. 8(b)). This type of pile-up was not detected by profile measurements through the centres of the neighbouring edges, where small sink-ins were typically observed.
The pile-up patterns described above were very pronounced for indentations close to the hard and soft orientations shown in Fig. 9 but more complicated behaviour was observed for intermediate angles. The additional contact depth due to pile-ups was less than 50 nm for indentation angles up to about 40°and then became more pronounced reaching a maximum of approximately 185 nm for orientations near 90°, as demonstrated in Fig. 8(b).
Pile-ups predicted by CPFEM simulations had a shape and position consistent with experimentally observed patterns although the maximum heights of the observed hills were smaller than in experiments and reached only 50 nm. This difference could be due to constraint from nearby crystals in the real sample enhancing the pile-up; when the dimensions of the simulated cube were halved the simulated pile-up profile was found to approximately double. Alternatively (or additionally) this discrepancy could be due to the approximations and simplifications made in the slip law. CPFEM was also able to predict sink-in behaviour for indentations away from the [0001] direction; however, the sink-in size observed in simulations was larger than the experimentally measured values.
Surprisingly, the observed patterns contradict the surface topography results of Tsuya (1967) from Vickers indentations of beryllium monocrystals: (0001) slip traces have been observed on the surface around indents in beryllium in a "soft" orientation. Surface traces perpendicular to (0001) were observed only at 300°C which the authors attributed to an increase in prismatic slip at this temperature. Detailed work by Zambaldi (Zambaldi et al., 2012) on the crystallographic dependence of pile-up locations in titanium also demonstrated that the slip systems with the lowest CRSS (< a > prismatic in the case of α-titanium) has the dominant influence on the indentation pileup behaviour. Our results demonstrate that while the CRSS of the prismatic slip is almost 5x higher than for basal slip, it determines the nanoindentation pile-up locations for soft beryllium orientations. Fig. 9(b) shows changes in contact area between the indenter and the sample due to pile-ups and sink-ins, where the contact area correction coefficient, C , is shown as a function of indentation angle. C is calculated as: The sink-in behaviour dominates for indents into grains with surface normals oriented close to the c-axis, reducing the contact area. The corresponding correction coefficient was measured be in the range of approximately −0.1 to −0.2. Pile-up dominates for indents into grains with surface normal orientated far away from [0001]. The correction coefficient gradually increases with indentation angle, reaching 0.25 for indents into grains with surface normal orientated perpendicular to [0001].
The averaged correction coefficient as a function of indentation angle between the c axis of the grain and the loading direction can be fitted to a linear function: Coefficients a 1 = 0.0035 and a 0 = −0.123 if the full data set is used (all beryllium grades, excluding indents close to grain boundaries), the coefficient of determination R 2 = 0.72. If the fitting = C f ( ) is done separately for the different beryllium grades, the less pure S-200-F and S-200-FH beryllium grades have a larger a 1 coefficient (up to 0.043) and slightly smaller a 0 coefficient (−0.181) implying deeper sink-in deformations for the c-axis indents and a more pronounced trend towards pile-ups in soft orientations.
The sensitivity of pile-up/sink-in behaviour in beryllium to crystallography and indenter orientation as captured by the correction coefficients, C, in beryllium is in contrast to the behaviour observed in isotropic materials, where the susceptibility to piling-up or sinking-in and its contribution to indentation contact area at a specific indentation depth can be approximated as a constant parameter as, for example, proposed for strain-hardened and annealed states of FCC copper (McElhaney et al., 1998). It is generally observed that a low (σ y /E) ratio and the absence of strain-hardening both promote the formation of pile-ups whereas the converse leads to the creation of sink-ins (Bolshakov and Pharr, 1998;Fischer-Cripps, 2004). Because of the complex three-dimensional profile of the strain and stress fields under the 3-sided pyramidal indenter tip, it is difficult to compare different orientations of the same materials using a single scalar value for E, σ y and the strain-hardening parameter. However the simulated stress field (Fig. 7) and GND density map (Fig. 6) penetrate deeper into the material in the hard orientation relative to the soft orientation consistent with a higher (σ y /E) ratio; and the plastic strain under the indent is also higher (Fig. 6) consistent with GND hardening.

Indentation-derived elastic modulus
The Young's modulus of industrial polycrystalline beryllium is around 303 GPa (Walsh, 2009). However, it is known that Young's modulus of beryllium is anisotropic and varies with crystallographic orientation from 341 GPa (parallel to [0001]) to 287 GPa (perpendicular to [0001]) (Walsh, 2009), as shown in Fig. 10 by the blue line.
Since the indentation modulus is obtained from a complex, three-dimensional, non-uniform stress state, the indentation modulus with indenter aligned along a given crystallographic direction usually differs significantly from the Young's modulus for uniaxial V. Kuksenko et al. International Journal of Plasticity 116 (2019) 62-80 deformation in that direction (Haušild et al., 2014;Vlassak and Nix, 1994). An anisotropic elastic material law was also used to perform finite element simulations and the predicted beryllium indentation moduli for indentation angles of 0°, 45°and 90°were 330 GPa, 320 GPa and 303 GPa respectively (green dashed line in Fig. 10). The analytical calculation of indentation moduli for monocrystals in the principal materials symmetry directions proposed by Delafargue and Ulm (2004), using anisotropic elastic constants from (Knezevic et al., 2013;Walsh, 2009) (2) and (3), show the reverse trend with indenter orientation, increasing from ∼275 GPa when = 0 to ∼340 GPa when =°90 (red data and line in Fig. 10). The measured pile-up and sink-in contribution coefficient C estimated using topography mapping profilometry was used to correct the indentation modulus E c using equation (14) derived from the combination of equations (2), (3) and (12): The corrected values (black points and line in Fig. 10) gives an increase in calculated modulus values in indents where sink-in dominates ( 0°) and lowers the modulus when pile-up dominates ( 90°). The corrected moduli do not have any noticeable dependence on indentation angle and chemical composition for the investigated samples and are scattered around the average value of 296 GPa. This is likely due to the assumption that only a single crystal was being probed during indentation, while during the experiment polycrystalline beryllium samples were tested. As demonstrated in Fig. 7, the field of elastic interaction between the indenter tip and beryllium samples is almost two orders of magnitude larger than the indentation depth and was significantly larger than the average grain size of the tested samples. This implies that the elastic response which determines contact stiffness from equation (3) for the effective modulus calculation is averaged over many crystallographic orientations and only weakly depends on the orientation of the indented grain. The measured average modulus of 296 GPa is in good agreement with the Young's modulus of industrial polycrystalline beryllium, which is between 292 and 304 GPa when measured in the hot pressed S-200-E grade (Christman and Feistmann, F.J., 1972).
It is important to emphasise that while our experimentally measured indentation modulus, E c , did not have a clear dependence on crystal orientation, the orientation did determine the localised plastic deformation and hence the contact area correction coefficient, C.

Hardness data correction
In this section, hardness values for all of the investigated grades will be extracted using the standard method (Oliver and Pharr, 1992) and accounting for the contact area correction. The measured pile-up and sink-in contribution coefficient C can be directly used for nanoindentation hardness correction (H c ) combing equations (1) and (12) V. Kuksenko et al. International Journal of Plasticity 116 (2019) 62-80 where H O P was extracted from the nanoindentation test data using the standard Oliver-Pharr method. However, only between 11 and 30 topography maps of indentation prints were made for each investigated grade, corresponding to 15-35% of indents at each grade. Two different contact area correction approaches were used. The first approach is based on the assumption that indents with similar angles should have a similar contact area correction coefficient. Equation (13) was used to estimate C, which was then used to determination the corrected hardness values H ( ) c using (15). This method is less accurate than the approach detailed above; however, it allows hardness data correction without requiring time-consuming experimental profilometry mapping of every indent.
The second approach was based on the elastic-modulus-based correction method which has been successfully applied in ferriticmartensitic steel (Heintze et al., 2016) and tungsten alloys (Beck et al., 2017). The method is based on the assumption, that as the contact area affects both hardness and indentation elastic modulus calculations (see equations (1) and (3)), then, in a material with known indentation modulus and Poisson's ratio, it should be possible to calculate the area correction coefficient from the measured modulus without having to image the indents at all. This data correction method can also be used for data from indents on or near to grain-boundaries, where no unambiguous crystallographic information is available. In this method, the area correction coefficient was deduced by combining equations (2), (3) and (12): where E O P is the indentation modulus obtained according to Oliver and Pharr (1992) without pile-up/sink-in considerations and E is the measured indentation elastic modulus of the material. As demonstrated in the previous section, the measured indentation elastic modulus of 296 GPa is in good agreement with the Young's modulus of industrial polycrystalline beryllium. This value was used for E in equation (16). The corrected hardness obtained with this method is denoted: Both C and C E were used for correcting the nanoindentation hardness data where the crystallographic orientation of the indented areas were confidently known (i.e. inside grains), and C E was used for data analysis of indents made close to grain boundaries. Fig. 11 compares hardness data from four beryllium grades before and after the corrections and shows that while H E exhibits higher scatter and suggest slightly smaller anisotropy relative to H C , both correction approaches showed reasonable agreement. The numeric values describing the observed trends are summarised in Table 1. Fig. 11 shows that the nanoindentation hardness as a function of indentation angle has a characteristic sigmoid shape with maximum and minimum hardness for indentation directions (test surface plane normals) parallel and perpendicular to [0001] respectively. Similar trends were observed in the micro-and macro-hardness tests performed in beryllium monocrystals by (Tsuya, 1967) and (Hill and Jones, 1961).
For ease of comparison, the experimentally observed relation between the indentation angle ( ) and the hardness data was fitted to a Boltzmann type sigmoidal curve: where hardness varies from H max to H min , 0 is the angle at which the hardness is halfway between H max and H min , and S describes the steepness of the curve. The fitting parameters together with the R 2 coefficient are shown in Table 1. Since a large number of indents had close to 90°, the fitting parameter H min is a good approximation of the average hardness for the soft orientation. Unfortunately, only limited data were available for indentation close to [0001] and given the observed scatter, the fitted H max parameter may not accurately represent the "hard" grain hardness. Fig. 11 and Table 1 demonstrate that high purity grades PF-60 and S-65 had almost identical hardness for the same crystallographic orientation, whereas S-200-F and particularly S-200-FH had noticeably higher hardness for soft orientations. In both high purity grades H min = 2.4 GPa while H min (S-200-F) = 2.8 GPa and H ( min S-200-FH) = 3.2 GPa. For indentation angles, < 50°, there was no explicit difference in hardness. Fig. 11 demonstrates that the difference in hardness distribution between PF-60 and S-65 grades ( Fig. 4(b)) is entirely due to the difference in sample textures. VHP and HIP-ed grades have a random texture, whereas in the rolled grade, grains are strongly aligned with their "hard" (0001) plane parallel to the samples surface. Table 1, shows our combined indentation/EBSD approach was very useful to detect the differences between beryllium grades. Averaging of nanoindentation data was not sensitive enough to reveal the difference and was affected be texture as in the case of PF-60 grade. The hardness difference should originate from impurity effects. The S-200-F has a higher content of oxygen, carbon, iron, aluminium, silicon and magnesium (see Table 1). Aldinger (Beryllium Science and Technology, 2012) demonstrated that the critical resolved shear stresses (CRSS) for both basal and prismatic slip in beryllium is proportional to c , 2/3 where c is the impurity concentration in wt%, for a wide range of impurity types. At the same time, purity does not have a significant effect on the flow stress of the pyramidal slip systems (Beryllium Science and Technology, 2012). The increased CRSS of < a > -type slip (while < c+a > is unchanged) could explain higher hardness and smaller anisotropy for the less pure grades, but detailed fine-scale chemical analysis of the S-200 grades is needed for confirmation. Fig. 12 shows the orientation dependence of experimentally-determined hardness values, as a function of position of the indention V. Kuksenko et al. International Journal of Plasticity 116 (2019) 62-80 Fig. 11. Nanoindentation hardness from 4 beryllium grades for different crystallographic orientation, calculated with the: (a) Oliver-Phar method (Oliver and Pharr, 1992), (b) corrected via profilometry data using equations (13) and (15)  Calculated using the standard Oliver-Pharr algorithm (Oliver and Pharr, 1992).
c Contact area corrected data using equation (15) which considers pile-up/sink-in contribution using the crystallographically dependent area correction coefficient calculated with equation (13). d Contact area corrected data using equation (17)  V. Kuksenko et al. International Journal of Plasticity 116 (2019) 62-80 direction within the unit stereographic triangle. While there was no distinct variation in hardness with angle between the [1010] and [1120] directions was found (in is in good agreement with our CPFEM results and the results of Tsuya (1967) on single crystals Vickers indentation), several localised minima and maxima were observed on the stereographic triangles (Fig. 12). This was mirrored by the very high scatter of results especially for < < 50 75 reaching ∼35% in the S-200-FH grade (Fig. 11). Tsuya (1967) did not report scattering of the Vikers indentation hardness data in single crystals. The fact that the scattering is most extreme in the grade with the smallest grain size could indicate influence of the neighbouring grains. On the other hand, the single crystal Knoop hardness data of Tsuya (1967) did show scatter reaching ∼25% when θ = 62°. Britton et al. (2010) also observed relatively high scattering (∼20%) of hardness around 55°with the c-axis in large grain α-titanium. This suggests the scattering observed in our work may originate (at least to some extend) from effects other than then influence of neighbouring grains. A detailed investigation of the heterogeneity of the dislocation and impurity distributions are required.
The H max /H min ratio obtained using the Oliver-Pharr approach is about 1.98 in the high-purity grades and 1.65 in the low purity grades, and after the area correction, they increase to 2.62 and 2.21 respectively. This is higher than in other common hcp metals. An H max /H min ratio between 1.25 and 1.5 was measured in titanium (Britton et al., 2010;Zambaldi et al., 2012) and α phase of a Ti-6Al-4V alloy (Viswanathan et al., 2005). H max /H min = 1.14 was observed in single-crystalline zinc (Vlassak and Nix, 1994). An H max /H min ratio value of about 1.17 was measured in single-crystalline magnesium (Somekawa and Schuh, 2016), but in this case the prismatic plane hardness was higher than the basal plane. The primary slip systems at room temperature in hcp metals are either < a > basal (Zn, Mg, Be) or first-order < a > prismatic slip (Ti, Zr) while the CRSS for < c+a > (pyramidal slip) is typically higher. However in beryllium the CRSS for < c+a > slip is significantly is higher than < a > (Beryllium Science and Technology., 2012) and this is reflected in the highest H max /H min ratio.
To our knowledge, this is the first nanoindentation study which compares different crystallographic orientation in beryllium, however some macro-and microhardness data obtained with 4-sided Vickers pyramid tip is available. The nanoindentation Berkovich tip was designed to have the same ratio of projected area to indentation depth as the Vickers indenter and both indenters have the same representative strain within the specimen material (Fischer-Cripps, 2004), therefore, the similar H max /H min ratios are expected from the two tip geometries (Rother et al., 1998). H max /H min ratios of 2.24 and 2.55 using Vickers indentation with 200g and 2 kg loads respectively were obtained by Hill and Jones (1961) from indentation of a single crystal beryllium of unspecified purity. A ratio of 2.78 was observed by Tsuya (1967) from Vickers hardness measurements of the hard (0001) plane and soft (1010) and (1120) planes in the low-purity single crystal beryllium (300g load). Both previous experiments with Vickers tips (Hill and Jones, 1961;Tsuya, 1967) are is in good agreement with our H max /H min ratio obtained from nanoindentation data, but only if the pile-up/sink-in contact area correction is performed. Considering this, it is important to emphasize that while the average nanoindentation hardness of the beryllium samples investigated here remained almost unchanged after the contact area correction, such corrections are important for description and comparison of nanoindentation hardness data between different grades or sample treatments.

Conclusion
Four industrial beryllium grades were investigated by nanoindentation and average hardness values of 4.7 GPa, 3.5 GPa, 3.9 GPa and 4.5 GPa were obtained for the PF-60, S-65, S-200-F and S-200FH grades respectively with a standard deviation of 20%. By combining nanoindentation, EBSD and CPFEM it was possible to separate out the influence of grade and crystal orientation on hardness. Grades with "soft orientation" close to [1120] or [1100] were most sensitive to the grade. High purity grades PF-60 and S-65 had almost identical hardness for the same crystallographic orientation with 6.4 GPa in "hard" grains (c-axis indentation) and 2.4 GPa in "soft" grains. Low-purity S-200-F and S-200-FH had similar hardness distribution for indentation angles less than 50°, and significantly higher hardness (17 %and 34% correspondently) for soft orientations. Hardness anisotropy was higher in the pure grades and ratio of H max to H min were 2.6 in the PF-60 and S-65 grades and 2.2 in low-purity grades.
CPFEM simulations demonstrated that GND distribution as well as elastic and plastic deformation areas highly depend on the indented grain crystallography. It is predicted that the plastic strain area is ∼7x larger than the tip penetration depth, giving a criteria for filtering quiai-single crystal hardness data. Our simulations demonstrated that the elastic interaction field between the indenter tip and beryllium samples is almost two orders of magnitude larger than the tip penetration depth and extend significantly beyond indented grains. This explains why no unambiguous link between the indentation modulus and the crystallography of the grain was observed in all the investigated grades.
Experiments and CPFEM simulations demonstrated that localised surface deformation around indentation prints (pile-up and sink-in) is highly crystallographically dependent: the sink-in behaviour dominated in hard grains (those indented close to the c-axis) while in soft grains (indented orthogonal to the c axis) pile-ups dominated. Topography investigations of the residual indentation prints demonstrated that when compared to the contact area calculated by the standard Oliver-Pharr method, the actual contact area is reduced by 10-20% in hard grains, whereas in soft grains the actual contact area is increased by approximately 25%. Evaluation of the hardness data and comparison with earlier studies demonstrated that the contact area correction is a key procedure for description of nanoindentation hardness between different grades or sample treatments of beryllium. Two different contact area correction methods were applied: an elastic modulus based correction and a profilometry derived correction. Both approaches showed good agreement with each other and with the data corrected by actual topography investigation of the residual indentation prints.
V. Kuksenko et al. International Journal of Plasticity 116 (2019) 62-80 which is also conveniently expressed as an × [8 3] matrix = B AJ 1 for a given Gauss point q.  This can be repeated for the second and third rows of F p to obtain the second and third columns of × F ( ) p respectively and hence T at that Gauss point. This process can then be repeated over all 8 Gauss points z q . There are two options which are worth noting. One option is to evaluate J , A and B at Gauss point coordinates (e.g. = ± z 1/ 3 i q ). However as the nodal values of F p are themselves Gauss point values this corresponds to obtaining T at the Gauss points of a virtual inner element whose nodes are at the Gauss points of the actual element. A final step is then required to extrapolate back to the Gauss points of the real element which can be achieved using  . The results are identical, provided the element is not distorted, but the second approach is faster as the last step of extrapolation is not required.