Crystal plasticity model for creep and relaxation deformation of OFP copper

ABSTRACT We demonstrate a dislocation density-based crystal plasticity (CP) model approach for simulating mesoscale deformation and damage. The existing CP framework is extended to be compatible with the oxygen-free phosphorous copper microstructure that is the focus of this study. The key aim is to introduce relevant plastic deformation mechanisms and to develop a failure model capable of depicting creep damage in the material. The effect of local variations in material is evaluated, and the model response is compared with experiments and characterisation. The basis of this work is CP material modelling, including grain orientation and size, obtained using electron backscatter diffraction and experimental test data of real relaxation test specimens. This will yield a realistic description of texture and grain shape and, ultimately, accurate stress–strain response at the microstructural level for further evaluation of performance with respect to material creep(−fatigue) damage.


Introduction
The concept of a final repository of spent nuclear fuel involves underground storage in capsules designed for long-term integrity in the storage environment [1,2].To protect the capsule against corrosion, it has a solid top layer (overpack) of oxygen-free phosphorus-doped (OFP) copper that is also used to withstand the in-service loads under residual heating by the contents.The resulting potential failure mechanisms include creep damage through nucleation and growth of creep cavitation and cracking [3,4].OFP copper has relatively good creep properties at modest temperatures, similar to the requirements of, e.g., the windings of electric motors and generators, and has also been previously studied for the intended application of nuclear waste canisters [4,5].In particular, the addition of about 50 ppm of phosphorus to oxygen-free high thermal conductivity (OFHC) copper improves both creep strength and creep ductility by strongly reducing the embrittling effect of sulphur at grain boundaries [6][7][8].
This paper describes a crystal plasticity (CP) model approach for simulating mesoscale deformation and damage.The main goal is to introduce a sufficient number of plastic deformation mechanisms and to develop a damage model capable of depicting creep cavity formation and propagation in copper material.A model based on dislocation density is utilised for the investigated doped copper material [9,10] with modifications.The damage model is incorporated into the modelling framework to cause degradation of the material intragrain and in the region of grain boundaries.

Experimental
The chemical composition of the currently investigated OFP copper was measured using optical emission spectroscopy technique and is presented in Table 1.
To demonstrate the propagating creep cavitation damage in OFP copper, samples from a specimen were creep tested under a uniaxial constant load of 95 MPa at 200°C up to 27,000 hours until failure.These samples were polished using a Struers Tegramin sample preparation unit, following standard metallographic practices.These samples were investigated with a Zeiss Ultra Plus scanning electron microscope (SEM), operating at 15 kV using an InLens signal for creep cavity images.For electron backscatter diffraction (EBSD) studies, samples from the gauge section were cut along three directions: two faces and one cross-section.The chosen sample was creep tested under uniaxial constant load at 300°C and 60 MPa for about 880 hours.These samples were mounted using conductive resin and ground and polished using a Struers Tegramin sample preparation unit following standard metallographic practices.To remove the final deformation layer from the sample surface, the sample was polished over 24 hours using a VibroMet® 2 vibratory polisher at 0.5 Hz using MasterMet 2 colloidal silica solution.Crystallographic orientations were measured using large-area EBSD.The EBSD was carried out using a Zeiss Gemini 2nd generation microscope, operated at 20 kV primary voltage and 15 nA probe current.The samples were tilted to 69.7° and kept at a working distance of 14-16 mm (depending on the height of the samples) to obtain the Kikuchi patterns.The EBSD patterns were captured using an 80 × 40 image resolution at a speed of 250 frames per second using an EDAX digiview camera.More than 2 × 2 mm 2 EBSD maps were recorded using a hexagonal grid and a step size of 1 μm.The collected EBSD data were then processed and analysed using TSL-OIM analysis v8.The maps were filtered using the confidence index (CI) standardisation and neighbour CI correlation (up to a single iteration).

CP modelling
A dislocation density CP model is formulated to describe the deformation behaviour of the investigated alloy.The deformation gradient is multiplicatively decomposed into elastic and inelastic contributions.
The inelastic contribution consists of plastic deformation by dislocation slip and damage part by crystalline damage systems.
where γ s and À N S are the slip rate and orientation tensor for a slip systems s.A total of 12 slip systems in {111}<110> slip family are included in the model.The damage model extension is based on strain-like formulation similar to references [10,11], where they are treated as cleavage planes.In the present case, however, due to the nature of ductile damage of copper-alloys, the crystalline {111}-planes are considered as damage planes rather than brittle-like cleavage planes, as was also suggested by Sabnis et al. [12].No additional shear mechanisms [12] are introduced as it is considered that shear accommodation of the damage planes opening on {111}-planes is accomplished by the slip systems.The choice of using damage as a strain like variable allows also to close the cracks mimicked with the damage model in the case that loading is reversed.This leads to a situation in which the material can partially regain its strength during compression, which can be considered as a smeared crack approach in general.Following this framework, damage rate _ δ d controls inelastic deformation, and the damage orientation tensor is defined as n s c � n s c opening of the damage planes, introducing anisotropic damage behaviour depending on orientation.
A viscoplastic Norton-type flow rule is used with isotropic and kinematic hardening contributions.
where τ s is the resolved shear stress, χ s is the backstress generating kinematic effect, τ s c is the isotropic shear resistance, K is the viscous parameter, and n is a strain rate exponent.The shear resistance in the basic form without damage consists of solid solution strengthening term τ 0 , a Hall-Petch contribution τ HP , and dislocation interaction and hardening term.The length of Burgers vector is denoted by b s .
The solid solution strengthening term τ 0 is so small in this case that it is set to zero and thus to be considered as included in the Hall-Petch contribution τ HP .The Hall-Petch term includes a temperature scaling with the use of effective shear modulus at different temperatures μ T ð Þ and a reference shear modulus at room temperature μð300KÞ.A constant Hall-Petch parameter K HP is used, and the effective grain size is denoted by d.
The nonlinear kinematic hardening follows the Armstrong-Frederick-type hardening model, with two model parameters C and D.
The interaction matrix a sj eff is not taken as constant but instead as an evolving one following the work of Monnet et al. [9].The interaction coefficients will, therefore, change with respect to dislocation density, which contributes to the saturation of the dislocation hardening effect at high dislocation density levels.This leads to a form: where ρ obs ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi P j�s ρ s dis q is the obstacle defect density and ρ ref is the reference dislocation density.Six FCC interaction matrix coefficients in a sj ref are used to distinguish between different types of dislocation [9].
The evolution of dislocation density is expressed as follows: where K mfp is a mean free path parameter, set as temperature dependent in the current model.Parameters K forest and K coplan define the average amount of obstacles passed by the moving dislocations before arresting.Dislocation annihilation is controlled by the annihilation distance G c , which is set constant in the present model.It can evolve as a function of strain rate and temperature [13,14] if chosen otherwise, e.g. by not setting K mfp as temperature dependent.
The damage rate is also expressed with viscoplastic flow rule so that whenever opening stress σ dc exceeds the damage resistance Y c , the damage process is initiated and begins to evolve.As the rate of damage can be negative when crack closure occurs, the damage strain is constrained to non-negative values (the crack/ void is virtually fully closed or crack/void has not formed in the first place).
where σ dc is the opening stress acting on the damage planes, Y c is the evolving damage resistance, K c and n c describe viscosity, and M is Mandel stress.The damage resistance is defined as a function of softening modulus H, degrading damage resistance when cumulative damage strain d c increases.It is assumed that slip localisation decreases the damage resistance and makes it more favorable nucleation site for creep cavities and subsequent cracks.Therefore, we also introduce coupling to cumulative plastic slip similar to references [10,11].
where σ c 0 is the initial damage resistance of undeformed and undamaged crystal and b c is the coupling parameter between slip and damage mechanisms.
Cumulative damage is defined asd In general, alternative models for damage/cavity growth and coalescence are available based on porous single crystal basis applied successfully to CP formulations [15].The current model with the damage-strain-based formulation, however, allows crack closure as mentioned earlier.
The coupling between slip and damage also influences the shear resistance of slip.Local cavities/cracks are expected to affect the slip resistance by softening mechanism to accelerate slip activity in the damaged zones.Softening to slip resistance is introduced by coupling with damage by modification on Equation 4. Following the approach descriped in [10].
Analyses were conducted with three different representative volume elements and one EBSD-based 2D mesh shown in Figure 1.The aim of the study was to build a model capable of capturing the deformation behaviour of the OFP copper, and for that, simplified RVEs are used.The representativeness of the model is tested with two simplified models comprising 64 and 512 orientations (grains).Simplified models are used to enable fast enough calibration and computation times.Creep and (creep) relaxation simulations are performed with a synthetic polycrystalline model containing 250 grains.A random grain orientation is assumed in the models, and the average grain size measured from the samples is used as the effective grain size in the models.Damage (cavity) formation in the models is analysed with SEM-EBSD-based 2D mesh.A similar EBSD-based approach was used with micromorphic strain gradient model in ref [16].
The 64-grain model (Figure 1a) was used to calibrate the model behaviour, while the 512-grain model and 250-grain polycrystalline model (Figure 1b and c) were used for the analysis cases.The EBSD-based model (Figure 1d) was used to demonstrate the model's ability to mimic the damage (cavity formation and growth) during the deformation.

Microstructural analysis
Representative microstructural image of the tested OFP copper material is shown in Figure 2. Material has dominantly random orientation, and no significant texture is observed.Grain size distribution of the material is relatively wide including noticeable amount of small and some large grains.The number of small grains is found to be highly influenced by the grain reduction caused by annealing twins.During the identification of CP model parameter, the twins are not explicitly included in the simulations.Therefore, the apparent grain size effect predicted by the model for grain sizes of 180 and 260 µm (in Figure 3) is almost negligible.However, when EBSD-based simulation domains are used, the annealing twins are a concrete part of the microstructure with the relevant misorientation from the parent grain.Deformation twinning is not included in the CP model as it is expected that applied strain does not significantly increase the number of twins.This aspect can be investigated in more detail in future by adding deformation twinning as a part of the model [17].

Mechanical response
The model's capability to present the mechanical tensile behaviour in several temperatures is presented and compared with the experimental observations of damage and failure in Figure 3.
Figure 3a shows the model's stress-strain response compared to experimental curves without damage at different temperatures.It is observed that the model agrees well with the experiments in the range of 25 −200°C.However, at higher temperature models, prediction of the stress-strain behaviour starts to deviate from experimental more noticeably.This is likely due to the intermediate temperature embrittlement (ITE) phenomenon or ductility trough behaviour that usually causes a significant drop of tensile elongation to failure and reduction of area at fracture in tensile loading at temperature range of 300-600°C [18].The model is presently not able to capture this phenomenon in a satisfactory manner.Figure 3d presents the behaviour of the EBSD-based 2D model.Due to the two dimensionality of the EBSD-based model hardening tends to be overestimated.This is handled by modifying the parametrization as shown in Table 2.The 2D EBSD-based models are used to compare the damage accumulation and behaviour with experimental characterisation, but the deformation behaviour is generally better modelled with 3D structure.
Figure 4 presents a comparison of the model prediction and experimental (creep) relaxation test.The model captures the experimental stress curve with satisfactory accuracy; however, it is observed that the first peak stress is higher in the model response by 10-15%.The difference can be partly explained with data acquiring rate (peak stress in the experiment might not be captured), but this needs to be verified in the future.In the simulated (creep) relaxation experiment, only constant strain rate is used due to the lack of deformation data with several strain rates, and thus, the calibration of kinematic effects is not fully taken into account at present.
Figure 5 shows the simulated and experimental creep test.The main differences between the model  and experiments are partly due to the vagueness in measuring the initial strain caused by the experimental set-up.Also, distinguishing the initial timeindependent (plastic deformation) strain and strain caused by time-dependent phenomena (creep) is not straigthforward and they are not separated in the model.Imitating this, however, is not the primary goal of the modelling approach.Instead, the strain rate evolution also differs for during the first few hours but most importantly the minimum strain rate is captured well with the model.Deviation of initial strain in model from experimental is partly due to the inaccuracies in the experimental set-up but also the accuracy of the model could be enhanced with more advanced characterisation/experimental techniques.Especially quantification and regularisation the slip    localisation and damage nucleation need more deatiled characterization.

Creep cavity formation and propagation
The diffusion-assisted and stress-directed flow of matter can result in local vacancy accumulation to cavity nuclei.In spherical nuclei, classical models assume stability at a critical diameter of 4γ s /σ, where γ s is surface energy at a grain boundary and σ is tensile stress normal to it [19][20][21][22].The theoretical work has indicated the minimum stable cavity size to be in the range of 2-5 nm, which is beyond the detection limit of usually applied electron microscopy [23,24].Easier  nucleation can be expected in cases of locally reduced surface energy, e.g. at boundaries to second phases like sulphides [19,25].In terms of the corresponding levels of normal stress, local increase would be predicted to facilitate easier cavity nucleation, but this stress could also be relaxed by the influence of the emerging cavities [19].
Regarding creep cavitation damage, there are two extremes in the roles of cavity nucleation and growth.One extreme is allowing continuous or initially continuous nucleation so that the growing number of cavities will significantly contribute to the overall total cavity volume.In the other extreme, all or nearly all nucleation sites are activated early, and the subsequent increase in the total cavity volume is mainly due to cavity growth (and later coalescence).Both types of cavitation have been observed in some materials [26][27][28][29].For example, in OFHC copper, partial segregation of sulphur at the grain boundaries provides several nucleation sites for cavities during creep exposure [29].Although, in OFP copper, the nucleationsupporting grain boundary sulphides have been reduced through phosphorus doping, resulting in suppressed nucleation and reduced number density of cavities, cavitation damage is still mainly located in grain boundaries [6][7][8]19].Figure 6 shows SEM images of creep cavities in the grain boundaries of OFP copper.

Damage accumulation in model
The damage evolution in the model and its comparison to characterisation are presented in Figure 7, showing the evolution of disorientation angle and damage in the crystal plasticity finite element (CPFE) model at different macroscopic strains (a: 10%, b: 12.5%, and c: 17.5%), where the disorientation angle signifies the magnitude of the rotation of the crystal orientation due to inelastic deformation.Figure 7d shows the EBSD plot near the damaged region and the contour plot of the disorientation angle from the mean grain orientation.This comparison indicates that the evolution of the disorientation angle is closely correlated to the damage in the CPFE model, signifying that at the location of damage, there is a large rotation of the crystallographic orientation.This is also verified by the experiments, in which the region near the damaged material shows a large difference from the mean grain orientation.
Figure 8b shows the normalised Taylor factor (TF) gradient magnitude [30,31].The TF [32] is estimated considering the loading direction in the CPFE model, mean grain orientation, and the slip systems.Here, the TF gradient shows the expected mismatch in the plastic deformation of the grain for a given loading direction.The TF gradient indicates possible sites of grain boundary damage, which can be correlated to the damage distribution in the CPFE result in Figure 8c plotted at 10% macroscopic strain.Due to the smaller computation effort in comparison to CPFE, the TF gradient analysis TFA model could possibly be used as a preliminary high-throughput method to evaluate susceptibility to failure.However, the magnitude of the damage cannot be correlated with the TF gradient.Furthermore, as observed in the CPFE damage distribution plotted at 12.5% strain (Figure 8d), the evolution of damage parameter inside grain cannot be predicted by the TF gradient analysis.

Conclusions
A CP framework was suggested for OFP copper material.The dislocation density-based formulation was used to describe dislocation slip deviated plasticity.A model extension was proposed to include crystalline-level damage that can cause material loss of strength and ultimately failure.The capability of the model to reproduce the experimental finding and predict the behaviour of the OPF copper was presented using a finite element solver.
The main findings in the characterisation, model development, and simulation are summarised below: • The material microstructural features, such as the grain size distribution and material orientation, as well as creep cavities, were characterised with high-resolution SEM-EBSD detector.In copper samples exposed to creep deformation, different levels of creep cavitation damage, including single cavities, chains of large cavities, and coalesced large cavities (microcracks), were discovered.The nucleation of damage primarily localises at the grain boundary zones, and the propagation of damage follows the grain boundaries.• A CP model that couples plasticity and damage was formulated.It was shown that the model is able to describe some of the most important damage phenomena in the material, such as grain boundary-type damage and shear localisation-related damage growth in the polycrystalline material.
• The model was used to investigate the effects of uniform grain size (the average grain size in the experiments) in the material definition of the 3D models, as well as realistic grain size with EBSDbased model, and heterogeneous strain hardening and strain localisation behaviour.In the future, the effect of grain size distribution and twin structures will be investigated also with 3D models in more detail.Furthermore, the incorporation of the effect of strain rate is also planned as the experimental data are foreseen to become available.

Figure 2 .
Figure 2. Representative SEM-EBSD image of the creep relaxation sample with orientation distribution figures and grain size distribution.

Figure 3 .
Figure 3. Models' stress-strain response (dashed lines) compared to experimental test curves (solid lines).180 µm curve is shown only in a).Curves produced without damage are shown in a), and curves with damage are in c), and average dislocation density evolution during the loading (in model) is presented in b).EBSD-based model stress-strain response is presented in d).

Figure 4 .
Figure 4. Polycrystalline model (creep) relaxation behaviour compared to experiential test at 50°C with a constant 0.5% strain.

Figure 5 .
Figure 5. (a and b) Polycrystalline model creep behaviour compared to experimental data.(c) Creep response of the copper material at 150°C with a constant load of 120 MPa.

Figure 6 .
Figure 6.(a) A faceted creep cavity in a grain boundary triple junction of OFP copper sample.(b) Several, potentially coalesced cavities partially covered by lid due to sample preparation in a grain boundary of OFP copper sample.

Figure 7 .
Figure 7. IPF-Z map of the CPFE model used in damage simulation.The evolution of crystallographic orientation (shown by disorientation angle) is compared with damage for (a) 10%, (b) 12.5%, and (c) 17.5% for the CPFE model.(d) The IPF-Z EBSD map of the damaged sample, which is compared with the variation in crystallographic orientation (shown by disorientation angle) at different locations with respect to the mean grain orientation.

Figure 8 .
Figure 8.(a) The original grain map (used in EBSD-based mesh).(b) The normalised magnitude of the TF gradient (in loading direction).(c and d) The damage resistance at 10% and 12.5% macroscopic strain, respectively.

Table 1 .
Chemical composition (wt ppm, bal.Cu) of the studied OFP copper material.

Table 2 .
List of parameters and initial conditions.