Characterisation of slip and twin activity using digital image correlation and crystal plasticity finite element simulation: Application to orthorhombic $\alpha$-uranium

Calibrating and verifying crystal plasticity material models is a significant challenge, particularly for materials with a number of potential slip and twin systems. Here we use digital image correlation on coarse-grained $\alpha$-uranium during tensile testing in conjunction with crystal plasticity finite element simulations. This approach allows us to determine the critical resolved shear stress, and hardening rate of the different slip and twin systems. The constitutive model is based on dislocation densities as state variables and the simulated geometry is constructed from electron backscatter diffraction images that provide shape, size and orientation of the grains, allowing a direct comparison between virtual and real experiments. An optimisation algorithm is used to find the model parameters that reproduce the evolution of the average strain in each grain as the load is increased. A tensile bar, containing four grains aligned with the load direction, is used to calibrate the model with eight unknown parameters. The approach is then independently validated by simulating the strain distribution in a second tensile bar. Different mechanisms for the hardening of the twin systems are evaluated. The latent hardening of the most active twin system turns out to be determined by coplanar twins and slip. The hardening rate of the most active slip system is lower than in fine-grained $\alpha$-uranium. The method developed in the present research can be applied to identify the critical resolved shear stress and hardening parameters of other coarse-grained materials.


Introduction
Hardening in metals is due to interactions within, and between, the slip and twin systems (Kalidindi, 1998). The nucleation and motion of dislocations and twin nucleation and migration occur when the critical resolved shear stress (CRSS) for these deformation mechanisms is reached (Devincre et al., 2008;Ojha et al., 2014). Since the seminal research of (Asaro, 1983), there has been significant effort dedicated towards the development of crystal based models for the inelastic deformation of engineering materials. Models of this type can contain a large number of material parameters, particularly if there are a number of different slip and twin systems that can be activated. A challenge then is to develop robust experimental procedures for the determination of these parameters. This either requires conducting a series of different experiments on polycrystalline materials covering a range of loading conditions and or stress states, or testing single crystals with different crystallographic configurations with respect to the loading direction (Franciosi et al., 1979). For a wide range of materials, particularly if they have limited ductility, each of these approaches might be impractical or even impossible to undertake. This is the case for the primary material of interest here -α-uranium (Inouye and Schaffhauser, 1969;Lander et al., 1994).
In this paper we demonstrate how the material parameters can be determined from a single uniaxial test on a coarse grained material by comparing the detailed strain and displacement fields within the individual grains determined using digital image correlation (DIC) with crystal plasticity finite element (CPFE) simulations. For α-uranium, this allows us to determine the 8 material parameters that describe the response of the different slip and twin systems that are active at room temperature. The approach can be readily generalised to other material systems.
Polycrystalline samples with columnar grains and arbitrary boundary conditions can be modelled using the crystal plasticity finite element (CPFE) method (Roters et al., 2018), which considers the anisotropic elastic behaviour, plastic deformation due to slip and twinning. Given the grain orientation, CPFE simulations provide the strain distribution in each grain (Dunne et al., 2007) and can be compared to digital image correlation (DIC) measurements (Lim et al., 2014).
The CPFE method has been widely used in conjunction with DIC measurements. For instance, DIC measurements can determine the local strain in polycrystalline samples, which is used for small scale CPFE simulations (Nellessen et al., 2015;Grilli et al., 2018). CPFE simulations are able to predict the strain localization (Irastorza-Landa et al., 2016) and lattice rotation (Irastorza-Landa et al., 2017b) in single crystals and oligocrystals, as detected by high resolution DIC (Guan et al., 2017). The DIC results can be correlated with the presence of geometrically necessary dislocations near grain boundaries (Roters et al., 2004). DIC has also been used to correlate strain concentrations with fracture paths (Khan and Marrow, 2009). As shown in the present study, the comparison between CPFE simulations and DIC measurements allows us to find the CRSS and hardening coefficients for the plastic deformation mechanisms in coarse-grained α-uranium.
In the CPFE method, the CRSS and hardening of the slip systems can be described by dislocationbased models, where the dislocation densities on the slip systems are the state variables (Grilli et al., 2015). An increase in dislocation density leads to hardening, usually described using a matrix detailing the strength of interactions between slip systems (Arsenlis and Parks, 2002). The dislocation density increase in α-uranium during deformation has not been systematically measured and few discrete dislocation dynamics studies appear in the literature (Behmer, 2018), which focus on the effect of twin thickness on dislocation pile-up and yield stress. The CRSS for slip and twinning have been identified for single crystal (Daniel et al., 1971) and fine-grained α-uranium (McCabe et al., 2010), obtained by hot or cold rolling. The interaction matrix between slip systems has been calibrated for fine-grained α-uranium (McCabe et al., 2010) using polycrys-tal simulations and the visco-plastic self-consistent (VPSC) homogenization scheme (Hutchinson and Hill, 1970;Lebensohn and Tomé, 1993). However, the CRSS and hardening have not been widely studied for coarse-grained α-uranium.
In the CPFE method, the interaction between twins and dislocations can be described using a CRSS for twin systems that depends on the accumulated shear on the slip systems (Abdolvand and Daymond, 2013a). The hardening rate of existing twins is also governed by the activity of coplanar and non-coplanar twins (Roters et al., 2010). This is based on the observation that, in FCC metals, deformation twinning initially produces only coplanar twins, while non-coplanar twin variants appear at larger strains (Kalidindi, 2001). This has led to the introduction of self hardening and latent hardening coefficients that couple the CRSS for twinning with the twin volume fraction (Abdolvand and Daymond, 2013a). These coefficients are set to zero in constitutive models describing fine-grained α-uranium, calibrated using stress-strain curves (McCabe et al., 2010) and twin volume fraction measurements (Knezevic et al., 2012). However, polycrystal simulations and fine-grained samples, in which thousands of grains are considered and the twin volume fraction is low (Inouye and Schaffhauser, 1969), are not ideal to study twin-slip and twin-twin interactions. The same holds true for neutron diffraction experiments (Earp et al., 2018;Grilli et al., 2019a,b), in which the measured lattice strain is averaged over a millimetre-sized gauge volume that contains many grains.
Coarse-grained α-uranium, obtained by casting, has a base centred orthorhombic crystal structure (Jacob and Warren, 1937;Lukesh, 1949;Kapoor et al., 2015) and grains reaching a size of several millimetres . It is ideal to study slip-twin interactions because of the propensity to form several twins after quenching or deformation (Cahn, 1951(Cahn, , 1953. This is because the CRSS for twinning is comparable to the CRSS for slip. Plastic deformation mechanisms in α-uranium are strongly influenced by temperature (Daniel et al., 1971;Lloyd and Barrett, 1966), strain rate (Huddart et al., 1980), grain size (Taplin, 1967;Inouye and Schaffhauser, 1969), grain shape (Taplin and Martin, 1963), loading conditions (Zecevic et al., 2016a), hydrogen content (Calhoun et al., 2015) and alloying elements (Yanzhi et al., 2015;Chen et al., 2017). At room temperature, the most active slip system is the so-called "wall" slip system, which has the lowest CRSS up to about 400°C, followed by the "floor" slip system (Zecevic et al., 2016b). Other slip systems exist, such as the "chimney" slip system, with two variants, and the "roof" slip system, with four variants (Calhoun et al., 2018;Calhoun, 2016;Calhoun et al., 2013), but are observed only at high temperatures (T > 873 K) when the temperature is close to the α → β phase transition (Daniel et al., 1971;Kapoor et al., 2015). The most commonly observed twin system is {130} 310 (Brown et al., 2009), but a small fraction of '{172}' 312 twins are also observed (Zhou et al., 2016). The '{176}' 512 twin system is observed only at high strain rate (Ho, 2012;Rollett et al., 1991).
As shown in the present study, tensile bars with a characteristic width of a few millimetres can be cut out from a larger polycrystalline plate made of coarse-grained α-uranium. The extruded area can be chosen such that only a small number of grains are present in the tensile bar aligned along the loading direction. High quality metallic surfaces, suitable for electron backscatter diffraction (EBSD) measurements, can be produced (Sutcliffe et al., 2018) and the grain orientation can be found. In-situ digital image correlation (DIC) measurements are carried out to measure the axial strain in each grain and the lateral displacement of the tensile bar. The stress-strain curve of 3 each grain is subsequently obtained and compared with the simulated curves obtained with the CPFE method. A nonlinear optimisation procedure, based on the Nelder-Mead algorithm (Nelder and Mead, 1965), is used to find the CRSS and hardening rates of different slip systems. The optimisation procedure can validate the values of 4 parameters against data from a tensile bar containing 4 grains. Two tensile bars are used to validate the model by comparing the measured and simulated strain field. Different models for the hardening of the twin systems are compared (Kalidindi, 2001).
The hardening of the most active slip system ("wall" slip) in the coarse-grained material turns out to be lower than for fine-grained α-uranium (McCabe et al., 2010). It will be shown that the hardening of the most active twin system {130} 310 due to coplanar twins and "floor" slip is able to explain the experimental results, and are therefore important mechanisms. This is possible by choosing a grain in which one variant of this twin system has a Schmid factor greater than for the slip systems. The experimental technique and computational method developed in the present research can be applied to identify the CRSS and hardening parameters of other coarse grained materials.
In Section 2 the crystal plasticity model for slip and twinning is described. Section 2.2 contains the three different models for hardening of the twin systems which were tested. Experimental details are reported in section 3. Section 4 contains the simulation results and the optimisation procedure to find the CRSS and hardening rates. The discussion and conclusions then follow in sections 5 and 6 respectively.

Crystal plasticity finite element modelling
The crystal plasticity framework is based on the decomposition of the deformation gradient F into an elastic part F e , describing the stretch of the crystal lattice, and a plastic part F p , describing slip and twinning (Roters et al., 2012): The plastic deformation gradient F p is an eigenstrain that does not contribute directly to the stress tensor. It has the following time evolution law (Kalidindi, 1998): The first term on the right hand side accounts for dislocation slip in the untwinned lattice, while the second term accounts for plastic deformation due to the creation of new twins. L p is the plastic velocity gradient which relates F p to its rate of change,Ḟ p .γ α (σ) is the dislocation slip rate on slip system α; defined by the slip direction, s α , and slip plane normal n α . f β is the twin volume fraction of twin system β (defined by the twin direction, s β , and twin plane normal n β ) which increases at a rate ofḟ β (σ) whereas the shear produced by the twin system, γ twin β , is a constant. Note, that bothγ α (σ) andḟ β (σ) are stress dependent, as discussed in section 2.1. The slip and twin systems used in the simulations are reported in table 1.
The slip directions and slip plane normals are rotated from the lattice reference frame to the sample reference frame by a rotation matrix R; which represents the grain orientation. R is updated at each time increment using the continuum elastic spin matrix, as explained in (Clausmeyer et al., 2011).
To model the nonlinear mechanical behaviour, at every time increment ∆t of the simulation, the Cauchy stress increment ∆σ = σ − σ 0 is calculated as a function of the deformation gradient, where σ 0 and σ are the Cauchy stress at the beginning and end of the time increment respectively. This is given by the anisotropic Hooke's law: where C is the fourth order elasticity tensor (Fisher and McSkimin, 1958) and W e is the continuum elastic spin (Belytschko et al., 2014). The elastic constants for α-uranium are reported in table 2. Details of the calculation of the elastic part of the small strain increment ∆ε e are reported in (Grilli et al., 2019a). To find the stress increment that satisfies equations (1), (2) and (3), an implicit equation for ∆σ is solved using a Newton-Raphson algorithm, as detailed in (Dunne et al., 2007;Dunne and Petrinic, 2006). The crystal plasticity framework is implemented in a user material subroutine (UMAT) for the finite element software Abaqus, which solves the equilibrium equation (Erinosho et al., 2013): ∇ · σ = 0.

Constitutive model for slip and twinning
In this section the constitutive model to calculate the shear strain rateγ α (σ) for each slip system and the rate of increase of the twin volume fraction for each twin system is reported. A  Table 2: Elastic constants (GPa) at room temperature for the orthorhombic structure of α-uranium (Fisher and Mc-Skimin, 1958;Beeler et al., 2013) in Voigt notation. power law relationship is used forγ α (σ) (Asaro and Needleman, 1985): whereγ 0 and n are constants that determine the strain rate dependence and rate sensitivity of the slip activity. τ α is the resolved shear stress and τ c α is the CRSS of slip system α, which depends on the dislocation densities. Similarly, the rate of increase of the twin volume fraction depends on the resolved shear stress τ β on the β twin system and is also assumed to follow a power law (Abdolvand and Daymond, 2013a): The CRSS for twinning, τ c β , is weakly temperature dependent (Christian and Mahajan, 1995) and it can depend on dislocation density (Abdolvand and Daymond, 2013a). The present model is not able to resolve discrete twins, but only the average volume fraction can be studied. τ c β has to be interpreted as the stress necessary to nucleate individual twins (Remy, 1978). The stress necessary for twin propagation (Qiao et al., 2016) and migration (Ojha et al., 2014) is typically lower than this value. These softening mechanisms, leading to stress drops (Patriarca et al., 2013), are not included in the present model.
A twin system can be activated only by a positive resolved shear stress (Kalidindi, 1998). The CRSS for slip is linked to the evolution of dislocation densities (Beyerlein and Tomé, 2008;Mecking and Kocks, 1981;Madec et al., 2002): where τ 0 α is the initial slip resistance, b α is the Burgers vector and µ α is the projected shear modulus (Grilli et al., 2019a) of the α slip system. ρ for α is the forest dislocation density, while ρ sub represents dislocation debris that forms stable dislocation substructures at the onset of stage IV hardening (Beyerlein and Tomé, 2008;Madec et al., 2002). The logarithmic term in (6) represents the stress for gliding dislocations to bow out through arrays of locked dislocations with spacing 1/ ρ sub . This model was originally developed by Beyerlein and Tomé (2008) for pure Zr, who provide further details.
The hardening rate is determined by the time evolution of the dislocation densities (McCabe where k 1 α are dimensionless parameters andd α is the annihilation distance for slip system α. The initial slip resistances τ 0 α for α = 1, 2, 3 ("wall", "floor", "chimney" slip systems) are three of the four parameters used in the optimisation procedure for the CRSS, described in section 4.1, and the dimensionless parameters k 1 α are three of the four parameters used in the optimisation procedure for the hardening.

Constitutive model for twin systems hardening
In the present study three hardening models for the twin systems are compared. The CRSS of the twin systems (130)[310] and (130)[310] will be indicated by τ c (130) and τ c (130) in the following. The corresponding twin volume fractions will be indicated by f (130) and f (130) . The first hardening model takes into account the interaction between twin and slip systems (Abdolvand and Daymond, 2013a). The CRSS for twinning is approximated with a linear function of the forest dislocation densities: The second hardening model (Roters et al., 2010) takes into account the interaction between the noncoplanar twin systems (130) The third hardening model (Roters et al., 2010) also includes the interaction between coplanar twin systems: Twin coplanar: τ 0 β is the initial twin resistance, µ β is the projected shear modulus (Grilli et al., 2019a) on the twin system and k 1 β is a dimensionless interaction parameter. τ 0 β is the fourth parameter used in the optimisation procedure for the CRSS, described in section 4.1. k 1 β is the fourth parameter used in the optimisation procedure for the hardening. The three hardening models are referred to as the twin-slip, twin-noncoplanar and twin-coplanar model from here on. The complete set of model parameters are reported in Table 3.

Experimental details
Experiments were performed to measure the full-field surface strains during tensile deformation of cast α-uranium specimens with a known microstructure, using digital image correlation (DIC). Two dog-bone specimens were manufactured by AWE plc via wire electro-discharge machining (EDM) from a cast plate. The gauge dimensions were 32.00 mm × 6.00 mm, with a specimen thickness of 1.50 mm.

Sample preparation and EBSD
Sample preparation and EBSD microstructural characterisation were performed at the Interface Analysis Centre (IAC), School of Physics, University of Bristol, UK. Each specimen was mounted in Struers ClaroCit resin, and mechanically ground with 120 grit SiC paper to remove the oxide layer. Grinding with progressively finer papers, from 320 to 2500 grit, was followed by polishing with 6 µm and 3 µm diamond suspension. The specimens were removed from the resin with acetone, and electropolished at a voltage of 12 V in an electrolyte comprising a 10:6:6 volumetric ratio of ethanol, ethylene glycol, and phosphoric acid. The electrolyte was continuously stirred and polishing proceeded for between 10 and 20 minutes, until a smooth, reflective surface was observed.
EBSD analysis was performed on a Zeiss Sigma HD VP Field Emission SEM, with EDAX EBSD camera and associated OIM software. A series of large-area maps covering the whole gauge area were recorded with a slight overlap between each scan region. The EBSD datasets were post-processed using the MTEX toolbox (Bachmann et al., 2011) in MATLAB. The individual orientation maps were stitched together to produce a single map for the gauge region of each specimen. Figures 1 and 2 show orientation maps from the two tensile bars, with colours representing the crystal direction parallel to the z direction (into-page). The microstructure is coarse, with sub-grains in excess of 200 µm in diameter clustered together into regions of similar orientations. Grain clusters are labelled according to their orientation; with 4 grains in the first tensile bar and 6 in the second.

In-situ mechanical testing
A speckle pattern of black and white acrylic paint was applied to the surface of each specimen with an airbrush. Tensile testing was performed on a Shimadzu Autograph AGS-X universal test machine with a 5 kN load cell, connected to a data acquisition PC running Shimadzu Trapezium-X software.
The specimen was enclosed in a bespoke sample holder to prevent contamination of the tensile grips. A window in the front face of the holder allows viewing of the specimen surface with the DIC cameras. The specimen was loaded in parallel with a thin ductile copper bar which would absorb energy in the event of specimen failure to prevent loss of radioactive material. The mechanical properties of the copper were characterised before the test; its flow curve was used to calculate the stress in the α-uranium specimen from the total load measured by the load cell. Each specimen was loaded in displacement control at a rate of either 0.05 mm/min or 0.025 mm/min. DIC was performed with a GOM ARAMIS 5M stereo DIC system. The camera mounting frame was rotated by 90 • from its usual horizontal orientation such that the affixed cameras were positioned above one another and angled inwards towards the specimen at an angle of 12.5 • . This allowed the gauge section of the specimen to be viewed in a portrait orientation. A calibration process allowed the co-ordinates in each image to be related to 3D co-ordinates in a small calibrated volume centred on the specimen surface. Two camera DIC systems as used here allow the separation of in-plane strains from out-of-plane displacements, which cause strain artifacts in single-camera DIC. Illumination was provided by polarised light from high-intensity LEDs. A pair of DIC images was automatically acquired in the ARAMIS software every 2.5 seconds during the test.
The DIC analysis was performed with a facet size of 20×20 pixels and a 50% overlap between each facet. A rigid-body movement of the specimen before testing was used to calculate the inplane strain error as ±0.1%, and the error in out-of-plane displacement as ±30 µm. The boundaries between the grain clusters in figure 1 were used to segment the DIC data such that the mean strain within each grain cluster can be calculated. This enables the plotting of stress-strain curves for the entire gauge volume (figure 4) and for the individual grain clusters (figure 5), for comparison with the simulation results.

Simulation details and results
The representative volume used for the first tensile bar is shown in figure 3 (a). The bar was 32 mm long and 6 mm wide. Displacement along the x axis was applied on the surface x = 32 mm. Zero displacement along the x direction was imposed on the surface x = 0. The point (x,y,z) = (0, 3, 0.75) mm was fixed to prevent translation of the surface x = 0. The points (x,y,z) = (0, 3, 1.5) mm and (x,y,z) = (0, 6, 0.75) mm can translate only along the z and y directions. This prevents rotation around the x axis of the surface x = 0. The texture of the representative volume consisted of four grains with different orientations ( figure 3). This is a simplification compared to the experimental EBSD map (figure 1) but is useful to maintain constant Schmid factors in each grain and, therefore, to better understand the slip-twin activity in each grain. This approximation also avoids the introduction of spurious data at points where the experimental EBSD has larger uncertainty. An uncertainty quantification on the grain orientations has been carried out and its effect on the computational results is discussed in section 5.
The coarser mesh used for the optimisation procedure is shown in figure 3 (b). It consisted of 434 elements with an average size of around 1 mm. The coarser mesh is necessary to reduce the computational time during the optimisation procedure, which requires hundreds of simulations. A finer mesh, of 19500 elements, shown in figure 3 (c), was used to repeat the simulation with the optimised parameters and compare the simulated and experimental strain fields.
The maximum time increment used in the simulations was 0.04 s and the total time was 20 s, after which a total strain of 2% was reached. Therefore, the simulated strain rate was 10 −3 s −1 : this value corresponds to the constantγ 0 in equation 4 and was chosen to simulate the regime in which the model is strain rate independent (Capolungo et al., 2009). In this regime, the resolved shear stress on a slip or twin system remains close to the CRSS τ c α or τ c β respectively. At the beginning of the simulations, the twin volume fractions were set to zero and the dislocation densities were set to 10 10 m −2 .

Optimisation procedure
The first optimisation procedure was carried out to find the initial CRSS parameters τ 0 α and τ 0 β for "wall" slip (α = 1), "floor" slip (α = 2), "chimney" slip (α = 3) and the {130} 310 twin (β = 1). A nonlinear optimisation procedure, based on the Nelder-Mead algorithm (Nelder and Mead, 1965), was used, which is implemented in the scipy.optimize package (Jones et al., 2001). This algorithm was used because it does not require an analytical expression for the Jacobian. At every iteration of the algorithm, the simulation with the coarse mesh in figure 3 (b) was run with a set of parameters τ 0 α , τ 0 β . For each grain in figure 1 the axial strain was averaged; denotedε G1 ,ε G2&3 andε G4 , for grains 1, 2 and 3, 4 respectively. The strain was averaged together in grains 2 and 3 because their mechanical behaviour was similar. The total strain in the tensile bar is referred to asε gauge and the stress was averaged over the load surface in figure 3 (a). This post-processing allows stress-strain curves to be extracted from the tensile bar and from the individual grains. The yield stresses of the tensile bar, σ gauge y , and of the individual grains, σ G1 y , σ G2&3 y , σ G4 y , were defined using the 0.1% offset criterion (Pham et al., 2013) for both the experiment and the simulation. The residual minimised by the Nelder-Mead algorithm was: where the subscripts exp and sim indicate experimental and simulated quantities. After every iteration of the Nelder-Mead algorithm, the parameters τ 0 α and τ 0 β are updated. This first optimisation procedure is carried out without considering the hardening of the twin systems because this does not significantly affect the yield stress. The CRSS parameters τ 0 α and τ 0 β that minimise the residual R CRSS in equation (13) are given in table 4.

CRSS parameters optimisation
Eqn. (6) and (9) Table 4: Parameters obtained from the optimisation procedures for the CRSS and for the hardening rate. The three values of the parameters k 1 1 , k 1 2 , k 1 3 , k 1 4 and k 1 β are obtained using the twin-slip, twin-noncoplanar and twin-coplanar model respectively.
A second optimisation procedure was then carried out to find the hardening parameters k 1 α and k 1 β for "wall" slip (α = 1), "floor" slip (α = 2), "chimney" slip (α = 3) and the {130} 310 twin (β = 1) using the CRSS parameters τ 0 α and τ 0 β obtained from the first optimisation procedure. The hardening parameters k 1 α and k 1 β were updated after every iteration of the Nelder-Mead optimisation algorithm, in which a simulation with the coarse mesh in figure 3 (b) up to 2% total strain was run. The maximum stress of the tensile bar reached during the deformation is denoted σ gauge . The maximum total axial strain in the tensile bar is denotedε gauge , and the maximum average axial strains reached during the deformation in the individual grains,ε G1 ,ε G2&3 ,ε G4 , were calculated after each simulation. The residual minimised during the second optimisation procedure: where and was determined for the three different hardening models for the twin systems reported in section 2.2. The hardening parameters k 1 α and k 1 β that minimise the residual R hard in (14) are reported in table 4.
The simulated stress-strain curves obtained using the optimised parameters and the three different hardening models for the twin systems are shown in figure 4 alongside the experimental data. Figure 5 shows stress as a function of the average axial strainsε G1 ,ε G2&3 ,ε G4 in the individ- ual grains for the three different hardening models for the twin systems. The stress-strain curves agree with the experiment for the twin-slip ( figure 5 (a)) and for the twin-coplanar (figure 5 (c)) interaction models. In grains 2 and 3, the maximum stress agrees with the experiment, as shown in figures 5 (a) and (c), while the maximum strain is underestimated by the models. The twin-noncoplanar interaction model is not able to match the stress value reached in the experiment, as shown in figures 4 (b) and 5 (b). It overestimates the strain accommodated in grain 4, as shown in figure 5 (b). This will be interpreted in the following section according to the activity of the slip and twin systems.

Fine mesh simulation results and model validation
Once the optimal parameters τ 0 α , τ 0 β , k 1 α , k 1 β were found, the simulation of the first tensile bar was repeated using the finer mesh in figure 3 (c). The model chosen for the hardening of the twin systems was the twin-coplanar model described by equation (12) as this provides the best optimised fit to the experimental data. The "wall" slip system was the most active because of the lower CRSS, as reported in table 4. Figure 6 shows the dislocation density ρ for 1 in the "wall" slip system at 1% total strain. "Wall" slip is mostly active in grain 1 and the dislocation density grows by a factor of approximately 100 compared with the initial value. Figure 7 shows the twin volume fraction f β at 2% total strain. It grows in the upper part of grain 4 and at the interface between grain 1 and grain 4, reaching a value of about 12%. We have Figure 6: Dislocation density ρ for 1 in the wall slip system in the first tensile bar at 1% total strain using the twin-coplanar model of equation (12). found that the contribution to the twin volume fraction of the system (130) 310 is much larger than the contribution of 1 30 31 0 . The twin-slip model and the twin-coplanar model lead to the correct prediction of the stress-strain curve in grain 4, as shown in figure 5 (a) and (c). The optimisation procedure allows us to find the interaction coefficients between coplanar twins and between slip and twin systems. The detailed strain and displacement fields within the grains provide additional information that can be used to validate the model resulting from the optimisation procedure. The axial strain concentrates in the upper part of grain 1 due to the activity of the "wall" slip system, as shown in figure 8. This agrees with the DIC measurement. The same slip activity leads to the concentration The measured and predicted out-of-plane displacement u z is shown as a surface plot in figure  10. The simulation results can be interpreted by studying the direction of the Burgers vectors of the "wall" slip system and of the {130} twin system. The largest component of these vectors is the u z component, which is positive. This is found by multiplying the rotation matrices R of grains 1 and 4 and the slip and twin directions s 0 α=1 and s 0 β=1 in the lattice reference system. Therefore, the large activity of the "wall" slip system in the upper part of grain 1 leads to positive u z in that region. The same can be stated for the large twin activity on the bottom part of the grain boundary between grain 1 and 4, shown in figure 7. This agrees with the DIC measurement in figure 10 (a).
The twin-noncoplanar model described by equations (10)-(11) does not predict strong hardening of the most active twin system (130)  An independent model validation was obtained by simulating the second tensile bar with 6 grains shown in figure 2. A fine mesh with 15984 elements was used. The boundary conditions were the same as for the first tensile bar, shown in figure 3 (a). The same model parameters τ 0 α , τ 0 β , k 1 α , k 1 β , optimised for the first tensile bar, were used. The model chosen for the hardening of the twin systems was again the twin-coplanar model described by equation (12). The measured and predicted axial strains are shown in figure 11. The "wall" slip system is active in grains 1 and 5, leading to strain concentration in those grains. The same slip activity leads to the concentration of the shear strain component ε xy in grain 1 and 5, as shown in figure 12.
The measured and predicted out-of-plane displacement u z is shown in figure 13. A rotation of about 1°around the x axis is added to the simulation results to match the experimental inclination on the surface x = 0, which was not present in the first tensile bar. Both the measurement and the simulation show a negative region of u z corresponding to grain 1.

Discussion
The optimisation method presented in section 4.1 shows that it is possible to use a single tensile bar experiment and DIC measurements to calibrate the yield strength and hardening rate of a crystal plasticity constitutive model for α-uranium. This material has slip and twin systems with very different properties. The reliability of the method is confirmed by the correct prediction of the strain and displacement fields, as shown in figures 8-13. This optimisation technique can be readily adapted to other materials. The stress-strain curves in figure 5 (a)-(c) show that the {130} twin system has a significant strain hardening response. For instance, the yield strength of grain 4 in the first tensile bar increases as the bar is loaded even though most plastic deformation in that grain is provided by a single twin variant (130) 310 , as shown in figures 6 and 7. Early studies on twin initiation and growth assumed a constant CRSS for twinning (Christian and Mahajan, 1995). Recent experimental studies (Patriarca et al., 2013) and molecular dynamics simulations (Ojha et al., 2014), however, show that both slip activity and twins affect the CRSS for twinning. This is consistent with the present work, in which the optimisation algorithm automatically identifies the extent of the twin-slip interaction and the interaction between coplanar twins. The comparison between the twin-slip model and twin-coplanar hardening models, reported in section 2.2, shows that both models can describe the mechanical behaviour of the first tensile bar. It is likely that the twin variant (130) 310 interacts with both slip systems and coplanar twins. The present approach is not able to distinguish between these two interactions, but it can provide values for the interaction coefficients k 1 β . The determination of the coefficients k 1 β was not possible in previous works on fine-grained α-uranium, based on tensile and compression tests of textured samples (McCabe et al., 2010), and on neutron diffraction (Brown et al., 2009). This is because these measurements include information from thousands of grains and are not able to discriminate between the slip-slip and slip-twin interactions. Therefore, the interaction coefficients between the slip systems and the {130} twin system were set to zero (see table 2 in (McCabe et al., 2010)).
Values for the twin hardening coefficients k 1 β of α-uranium are not available in the literature, but a quantitative comparison with other materials can be made. (Abdolvand and Daymond, 2013b) used a Voce hardening law to describe Zircaloy, in which the CRSS for tensile twins was given in the small strain approximation: where Γ is the accumulated shear on all slip/twin systems and θ 0 β = 50 MPa. Assuming that the {130} twin system is the only plastic deformation mechanism, equation (12) for the twin-coplanar hardening model becomes: where the stress is expressed in MPa. This shows that the values obtained for k 1 β have the same order of magnitude as in other materials. Generally, the hardening rate found in the present experiments is lower than that found for fine-grained α-uranium. For instance, the dislocation multiplication prefactor of the "wall" slip system reported by McCabe et al. (McCabe et al., 2010) is k 1 1 /b 1 = 200 µm −1 while the present optimisation algorithm finds a maximum value of k 1 1 /b 1 = 43 18 µm −1 , as reported in tables 3 and 4. This is consistent with the finding that the average dislocation density grows faster in fine-grained materials (Haouala et al., 2018;De Sansal et al., 2010).
Mesoscale models for FeCr single crystals, based on molecular dynamics simulations (Ojha et al., 2014), show that the twin migration stress grows linearly with the residual Burgers vector left at the twin interface during twin-dislocation interaction (see figure 14 in (Ojha et al., 2014)). Those simulations are carried out in a single crystal with a characteristic dimension of 500w, where w is the dislocation core width. A migration stress of 167 MPa is found if the residual Burgers vector is |b r | = a, where a = 0.286 nm is the lattice constant (Ojha, 2013). Assuming w ≈ a, the residual Burgers vector in those simulations corresponds to a dislocation density ρ ≈ 1/ (500a) 2 ≈ 49 µm −2 . The linear relationship between the dislocation density and the twin migration stress gives an increase of 167 MPa / 49 µm −2 ≈ 3.4 MPa / µm −2 . Also equation (9) predicts an increase of the CRSS τ c β , above which twin volume fraction grows, that depends linearly on the dislocation density; the factor of proportionality is k 1 β µ β b β b α ≈ 14 MPa / µm −2 . This comparison shows the consistency between the present model, calibrated without imaging individual twins, and atomistic models, which determine the interaction between dislocations and twins at the sub-micron length scale.
The uncertainty in the present optimisation procedure arises from the approximation that the four grains of the first tensile bar in figure 1 each have uniform crystal orientations. An uncertainty quantification procedure was carried out: the Euler angles representing the crystal orientation were sampled at different points in the 4 grains; for instance, this leads to an uncertainty of about 20% on the Schmid factor of the "wall" slip system. The optimisation procedure has been repeated with different Euler angles. The uncertainty is small, about 2%, on the constant friction stress τ 0 1 and on the dislocation multiplication prefactor k 1 1 of the "wall" slip system, and it is about 30% for τ 0 β and k 1 β of the twin system. The constant friction stress of the "chimney" slip system has a large uncertainty: if the value is changed, the residual R CRSS of the optimisation procedure does not change significantly. This is because of the low Schmid factor of the "chimney" slip system in the 4 grains and the higher CRSS.

Conclusions
An optimisation procedure was developed to find parameters of a dislocation-based constitutive model, implemented in a crystal plasticity finite element framework. DIC measurements were made during tensile tests on coarse-grained α-uranium bars containing a small number of grains. This allows us to find the average strain in a grain and to plot the stress-strain curves of individual grains, which can be directly compared to the simulations.
The constitutive model includes the activity of 8 slip systems and 2 twin variants, as previously identified in α-uranium (Zhou et al., 2016). Three different models for the hardening of the twin systems were compared. Models which account for the interaction between twin and slip systems and the interaction between coplanar twins can explain the experimental data. The optimisation procedure is able to find the values of the interaction coefficients, which have not been identified in previous studies on α-uranium.
The optimisation procedure is based on consideration of the average response of the grains within the coarse-grained specimen. The resulting model has been validated by comparing the detailed experimental and simulated strain and displacement fields within the same specimen and for a second tensile specimen with a different grain structure. These are interpreted by analysing the activity of twin/slip systems and the direction of their Burgers vectors.
The method presented can be applied to find the CRSS and hardening properties of plastic deformation mechanisms of other coarse-grained materials. The method has the advantage that the preparation of single crystal samples is not required and a large number of parameters can be determined from a single test.