A Continuum Dislocation Dynamics Crystal Plasticity Approach to Irradiated Body-Centered Cubic α-Iron

Ferritic-martensitic irons are used in the reactor pressure vessels (RPVs) of existing light water reactor nuclear power plants, due to their resistance to irradiation damage [1]. These materials are also under consideration for some advanced small modular reactor designs [2]. Despite this resistance, over extended exposure to radiation doses, these body-centered cubic (BCC) irons demonstrate irradiation damage, the primary radiation damage defect being self-interstitial atom (SIA) loops [3]. These irradiation defects interact with dislocations, increasing the brittleness of the material by impeding the movement of dislocations [4,5]. The impeded dislocation motion produces increased yield stress and significant hardening in the plastic regime, causing the material to behave in a brittle manner [2]. In addition to significant radiation loads, the RPV is subjected to thermal and pressure loading from power plant operation. The combination of these effects causes the degradation of the RPV over time as impeded dislocation movement in the microstructure reduces the ability of the material to accommodate mechanical loading [6]. RPVs must have sufficient toughness to avoid fracture during normal operation or transient events, and embrittlement reduces material toughness over time.

Crystal plasticity models play a key role in capturing how mesoscale changes in the microscale affect the engineering scale material properties. Recent developments in crystal plasticity modeling have pushed toward models based on dislocation density [7,8]. These models, developed for face-centered cubic (FCC) materials, use dislocation density evolution equations with origins in the simple generation and annihilation terms proposed by Kocks et al. [9]. These models separate out populations of screw and edge dislocations; yet, the evolution equations for the dislocation populations still reduce to a binary multiplication and annihilation balance. Also working with FCC materials, Roters et al. and Ma et al. used the concept of mobile and immobile dislocations in their dislocation evolution models [10,11].

Radiation damage has been incorporated into FCC crystal plasticity models through dispersed barrier models [12,13]; in these models, the radiation damage defect density is held constant. Newer crystal plasticity models focused on RPV applications couple irradiation defects to the dislocation evolution model by adding to the slip resistance calculation [14–17]. Among these crystal plasticity models that couple dislocations and SIA loop evolution, two employ a dislocation evolution model that tracks both mobile and immobile dislocations. These dislocation evolution models include terms for observable dislocation interaction mechanisms, including locking and cross slip [15,17]. The implemented cross slip term, however, relies on an approach that can saturate the dislocation density and does not account for the stochastic nature of cross slip within the physical crystal microstructure.

Several groups have focused on the problem of modeling the interactions among dislocations and SIA loops in BCC metals; yet, none have fully captured the physical interactions among key microstructural components, particularly in regard to cross slip dislocations. Since this work focuses on capturing the physical interactions among dislocations and SIA loops, we performed simulations using α-iron exclusively.

In this work, we present a dislocation-density-based continuum crystal plasticity model composed of physically based interaction mechanisms in the dislocation evolution equations, including a stochastic approach to cross slip calculations in order to capture the random nature of this dislocation mechanism. The dislocation evolution model is coupled to an equation for SIA loop evolution, and hardening is calculated as a function of both dislocation and SIA interactions. Single-crystal verification experiments are used to calibrate the continuum dislocation dynamics (CDD) crystal plasticity model prior to assessing the model’s predictive capability using a simplified polycrystalline geometry. In Sec. 2, we introduce the constitutive model as well as details on the dislocation and SIA loop evolution equations. In Sec. 3, we describe our calibration of the CDD crystal plasticity model to single-crystal α-iron tensile experiments. In Sec. 4, we present the results of the finite element simulations, including a discussion of cross slip approaches and a comparison to lower length scale approaches.

The challenge in continuum level models is to adapt the discrete probability model for continuum dislocation density values while retaining the physical foundation of the cross slip model. This work explores the effects of two different models for stochastic cross-slip representation within a continuum framework: a stochastic-only model and a stochastic Monte Carlo combination model for calculating dislocation cross slip. The stochastic-only model is taken from Patra and McDowell [15] and denoted in this work as “stochastic-PM” the combination stochastic Monte Carlo method introduced in this work is given the notation “stochastic-MC.”

Within the stochastic-PM approach, Eq. (19), all slip systems participate in each cross slip transfer of dislocations: each slip system loses dislocations and gains cross slip dislocations from all other slip systems in the cross slip family at every simulation time-step. In contrast, in the stochastic-MC approach, Eq. (20), all of the cross slip dislocations from the first slip system are assigned to a single Monte Carlo selected slip system, all of the dislocations from the second slip system are assigned to a single selected receiving slip system, and so on. This procedure is performed at all individual material integration points at the simulation time-step.

This crystal plasticity model has been implemented in the Multiphysics Object Oriented Software Environment (MOOSE), which was developed at Idaho National Laboratory to solve coupled physics simulations in a finite element method framework [33]. We calibrated our proposed crystal plasticity model by focusing first on unirradiated α-iron single-crystal data before verifying the implementation of the irradiation-defect-hardening models against lower length scale simulations.

Using single-crystal tension test experimental data [34], we calibrated the coefficient parameters for the six mechanism interaction terms used in the dislocation evolution models, Eqs. (9) and (10), against α-iron single-crystal data from three different loading directions [35]. These particular loading orientations change the influence of the cross slip term (α₅ in Eq. (15)) by varying the number of active slip systems: [001] has four, $[01¯1¯]$ has two, and $[3¯48]$ has only a single favored active slip system.

The single-crystal calibration simulations were performed on a 1 mm³ cube mesh consisting of 216 8-node hexahedron elements. A mesh convergence study was performed, in all three loading directions shown in Fig. 1, with a series of 8, 64, 216, and 512 elements. The loading direction with the greatest number of active slip systems, [100], demonstrated no mesh dependence while the $[3¯48]$ loading direction demonstrated the largest mesh dependence. The $[01¯1¯]$ and $[3¯48]$ loading directions demonstrate acceptable mesh convergence with 216 elements. Symmetric boundary conditions were used, and a displacement loading rate corresponding to a strain rate of 3.3 × 10⁻⁴ 1/s was applied, matching Keh [34]. The values of the elastic properties, glide velocity, and dislocation cross slip were held constant throughout this calibration process and are listed in Table 1. We assumed equal values for both the initial mobile dislocation density and the initial immobile dislocation density, set at 2.5 × 10⁵ mm⁻² [35] and distributed evenly across all slip systems.

The values of the dislocation mechanism terms, α₁–α₆, were fit using insights deduced from dislocation dynamics simulations [28]. Based on the results of these dislocation dynamics simulations, we gave priority in the fitting to the first two terms in the mobile dislocation evolution rate, Eqs. (11) and (12), and adjusted the values of the coefficients to obtain agreement with the experimental curves for each loading orientation. The crystal plasticity simulations using the finalized parameter values are presented in Fig. 1 and compared against the experimental data from Keh [34]. The specific values of the dislocation evolution parameters are listed in Table 2. The different Peierls stress values for the different loading directions was also observed by Taheri and Zbib [29]. The variance of the cross slip evolution parameter indicates that the role of cross slip is more significant in the deformation of the single active slip system loading direction compared to the loading direction with multiple active slip systems, as shown previously by Li et al. [28].

Inclusion of a cross slip term in the dislocation evolution equations is key for capturing the stress-strain behavior of the single slip loading orientation. Cross slip of dislocations away from the activated slip system mitigates the growth of this slip systems’ dislocation density; thus, cross slip prevents over-hardening of the effective stress response.

A key feature of our crystal plasticity model is the inclusion of the combination stochastic-MC cross slip term, Eq. (20), to help capture the random nature of the physical cross slip dislocation movement. Both the stochastic-MC approach and the more commonly applied stochastic-only-PM approach use the same calculation for the probability of cross slip, Eq. (18), which is a function of the applied shear stress. The difference between the two approaches lies in the method of determining which slip systems will interact during cross slip. Based on previous studies [28], we anticipated that the largest difference in the cross slip approaches would occur in the $[3¯48]$ loading direction. In our comparison of single-crystal simulations, we applied the same symmetric boundary conditions and loading rate, with the material parameters given in Tables 1 and 2, for the $[3¯48]$ loading direction. In the stochastic-PM cross slip approach, we adjust the leading coefficient, α₅, to produce parity in the density of the cross slip dislocations in the two approaches at the beginning of the simulations. With these consistent trends in the probability of cross slip among slip systems, Fig. 2, we examined the impact of the two cross slip approaches on the mobile dislocation density, Fig. 3. In both approaches, Figs. 3(a) and 3(b), the cross slip of dislocations away from the active $(2¯11)[111]$ system—as well as, to a lesser extent, away from the secondary activated $(1¯21)[111¯]$ and $(1¯21)[111¯]$ systems—relieves the growth of dislocations on the active slip systems.

The stochastic-PM approach, Fig. 4(a), introduces an apparent saturation limit in the $(2¯11)[111]$ system mobile dislocation density. The saturation from the deterministic approach reduces the mobile dislocation density of the $(2¯11)[111]$ system to such an extent that the mobile dislocation density of the secondary $(1¯21)[111¯]$ system grows to parity with that of the primary $(2¯11)[111]$ system. This saturation of the $(2¯11)[111]$ system’s mobile dislocation density runs counter to the expected continued dislocation density growth on the primary activated slip system. While the stochastic-PM approach does produce a stress-strain curve similar to the experimentally measured curve, as indicated by Fig. 4(a), the inability of this approach to maintain the $(2¯11)[111]$ system as the primary activated slip system demonstrates the unsuitability of the stochastic-PM approach to a physically based crystal plasticity model.

In contrast, in the stochastic-MC approach, Fig. 4(b), this transfer of dislocations via cross slip enables the dislocation density on the active $(2¯11)[111]$ system to grow at a moderate pace. Even though the mobile dislocation density on additional slip systems, $(1¯21)[111¯]$ and $(11¯2)[1¯11]$⁠, increases, the favored slip system, $(2¯11)[111]$⁠, consistently maintains the highest dislocation density value, as is expected in a single slip loading orientation.

Inclusion of anisotropy is necessary to capture the single slip system behavior of the $[3¯48]$ loading direction orientation. Similar cross slip probabilities are calculated on both the $(2¯11)[111]$ system and the $(1¯01)[111]$ system, yet the anisotropy of the slip systems, Eq. (7), allows mobile dislocation growth only on the softer {112} system. Since these two systems have the highest probability of receiving cross slip dislocations within the [111] cross slip family, it is therefore likely that dislocations will cross from $(2¯11)[111]$ to $(1¯01)[111]$ and $(1¯10)[111]$⁠, as shown in Fig. 5(a). Transfer of dislocations from the active $(2¯11)[111]$ system to the $(1¯01)[111]$ system relieves the dislocation growth on the active system; this dislocation growth mitigation prevents over-hardening of the crystal stress response. The interaction between the anisotropy and the stochastic-MC cross slip approach captures the single active slip system behavior expected under $[3¯48]$ loading.

The higher slip system resistance allows the $(1¯01)[111]$ and $(1¯10)[111]$ systems to absorb the cross slip dislocations from the $(2¯11)[111]$ system. This absorption of cross slip dislocations from an active system by an inactive system enables our crystal plasticity model to capture the near-ideal plastic behavior of the single crystal loaded in the $[3¯48]$ direction.

Verification of the SIA loop defect model consists of comparing the stress-strain curves generated by our crystal plasticity model to those from dislocation dynamics simulations [16,37]. As with the dislocation dynamics simulations, these crystal plasticity simulations were performed on 1 μm³ cubes. These cubes were loaded in tension in the [100] direction, with traction-free boundary conditions on the lateral sides. In alignment with the dislocation dynamics simulations, a displacement loading rate equivalent to a strain rate of 100 1/s was applied to the top of the cubes. Six different verification simulations were run, each with a different initial SIA loop density: 1.63 × 10¹³ mm⁻³, 8.15 × 10¹² mm⁻³, 3.61 × 10¹² mm⁻³, 1.63 × 10¹² mm⁻³, 8.15 × 10¹¹ mm⁻³, and none (unirradiated case). As per Barton et al. [16], an initial dislocation density of 2 × 10⁻⁷ mm⁻² was assumed (split evenly between mobile and immobile dislocations), and a Peierls stress of 80 MPa was applied. The other SIA-loop-specific parameters are given in Table 3. The remainder of the model constants used in these verification simulations is given in Table 1, and the dislocation evolution parameters for the [100] loading direction are listed in Table 2.

The results of these simulations compare reasonably well with published dislocation dynamics simulations [37] as shown in Fig. 6. As in Barton et al. [16], the use of anisotropic elasticity in the crystal plasticity model requires comparing the plastic strains in the dislocation dynamics simulations with those in our crystal plasticity simulations. Our model relies on a scalar form of SIA loop evolution, while Barton et al. employ a tensorial form of the SIA density rate to account for the 3D nature of SIA loop interaction [16]. The 3D interaction among the SIA loops and the dislocations is accounted for in our model with a cube root term, Eq. (21). Comparison of this crystal plasticity model with the dislocation dynamics simulation results [16] (Fig. 6) demonstrates alignment of the trends. The misalignment between the two model predictions is due in part to difference in the discrete and continuum approaches of the dislocation dynamics and crystal plasticity frameworks, respectively, with the continuum crystal plasticity method demonstrating less variability than the discrete dislocation dynamics model results. Based on these results, our scalar SIA loop density evolution rate model, Eq. (23), can acceptably replicate the lower length scale trends for varying initial SIA loop densities in a less computationally intensive manner.

We compare the results of our CDD crystal plasticity model to irradiated polycrystalline experimental data in order to evaluate the model’s ability to predict RPV behavior after exposure to radiation. As a first step, we compared a polycrystalline application of the CDD crystal plasticity model to experimental data on unirradiated α-iron. Our simplified polycrystalline geometry consisted of 27 equally sized cubic grains with a diameter of 250 × 10⁻⁶ mm [39]. Each grain was meshed with 216 elements, for a total of 5832 8-node hexahedron elements in the full mesh. The orientations of the grains were determined through random assignment of the three Bunge Euler angles, within the usual angle bounds, using the python random number generator with a normal distribution [40].

The dislocation evolution parameters from Table 2 for the [100] loading orientation were used in these polycrystalline simulations. The same α-iron material parameters from Table 1, except for the Peierls stress value: here we used the polycrystalline value of 11 MPa [23]. We applied an initial dislocation density value of 5 × 10⁷ mm⁻², which falls within the experimentally measured 7 ± 2 × 10⁷ mm⁻² range [39]; this value of initial dislocation density was selected by calibrating to the unirradiated data for polycrystalline α-iron. The initial dislocation density value was split evenly between the mobile and immobile initial dislocation densities. We applied the symmetry boundary condition to the model, with a displacement loading rate corresponding to the strain rate of 2 × 10⁻⁴ 1/s in order to match the experimental loading conditions.

The same mesh, Bunge Euler angles, and parameters were retained for the irradiated polycrystalline simulation. To capture the effect of the irradiation defects, we included the terms for the SIA loop evolution, Eq. (23), and interaction with the dislocations, Eq. (21). The parameters for these equations, selected to correspond to an irradiation dose of 0.1 dpa, are given in Table 4. We set the value of the dislocation slip-system hardening coefficient, α_SIA, to 0.6, which lies between the values of 0.37 and 0.7 used in other studies of irradiated α-iron [17,38]. Following the approaches of Li et al. and Chakraborty and Biner, we also lowered the value of the SIA loop annihilation coefficient by a factor of 100 [17,28]; we hypothesize that this reduction of the SIA loop annihilation constant is a result of greatly reducing the loop diameter from that used in the comparison with the dislocation dynamics simulations in Sec. 4.3.

The stress response of the cubic polycrstalline simulation is calculated as an effective von Mises-type stress measure, equally averaged across all quadrature points in the mesh. In Fig. 7, this effective second Piola–Kirchhoff stress measure is plotted against an effective strain, then compared to the experimental data. The CDD crystal plasticity simulations agree with the unirradiated experiment and irradiated experimental data trends. The CDD simulation results’ slight underprediction of the experimentally measured data in Fig. 7 indicates that the CDD crystal plasticity model lacks a hardening contribution. A possible correction to the underpredicted hardening is a strain-gradient-type term; such a term is often used in crystal plasticity frameworks to capture the effects of grain boundaries on the stress response of the bulk crystal [22].

The general alignment of our CDD model with the measured polycrystalline response—even when applied to a simplistic cubic geometry—demonstrates the capability of our crystal plasticity model to predict the hardening behavior of α-iron exposed to irradiation. The physically based dislocation and dislocation-SIA loop interaction terms of the CDD crystal plasticity model can successfully capture the complex physical mechanisms of irradiated α-iron.

We have developed and applied a CDD crystal plasticity model with dislocation evolution terms based in physical dislocation interaction mechanisms. The dislocation evolution is coupled with the SIA loop evolution in acknowledgement of the interstitial loops’ significant impact on the irradiation behavior of RPV steels. This model leverages the results of lower length scale molecular dynamics and dislocation dynamics simulations to establish evolution equations for both the dislocations and SIA loops. We calibrated the dislocation evolution components of the crystal plasticity model against single-crystal α-iron tensile experiments, with insights from dislocation dynamics studies.

The combined stochastic-Monte Carlo dislocation cross slip model, in conjunction with the anisotropic strength of the BCC slip systems, is necessary to correctly capture the dislocation behavior in a loading orientation selected for single slip system activation. The anisotropy correction factor on the {112} type slip systems allows these slip systems to absorb dislocations, without increasing their own mobile dislocation densities, from active lower length {111} type slip systems. The transfer of dislocations through cross slip from softer to harder slip systems prevents overhardening of the single-crystal response.

Verification of the coupling between the dislocation evolution and the SIA loop evolution was performed by comparing the simulation trends from this model with trends from dislocation dynamics simulations. We then applied the CDD crystal plasticity model to a polycrystalline application with a simplistic geometry. The CDD model demonstrates the capability to predict the stress response of irradiated polycrystalline α-iron, and the results of this mechanism-based CDD crystal plasticity model can be used to inform engineering scale models that rely on the stress response from an evolving microstructure under radiation and deformation conditions.

This manuscript includes work performed by some of the authors while contractors of the U.S. Government under Contract DE-AC07-05ID14517. This research made use of the resources of the High Performance Computing Center at Idaho National Laboratory.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.