## Abstract

Radiation-induced embrittlement of reactor pressure vessel (RPV) steels can potentially limit the operating life of nuclear power plants. Over extended exposure to radiation doses, these body-centered cubic (BCC) irons demonstrate irradiation damage. Here, we present a continuum dislocation density (CDD) crystal plasticity model to capture the interaction among dislocations and self-interstitial atom (SIA) loops in α-iron. We demonstrate the importance of modeling cross slip using a combined stochastic Monte Carlo approach and the role of slip system strength anisotropy in capturing stochastic cross slip interactions. Through these captured interactions, the CDD crystal plasticity model can capture both the stress response and the physical evolution of dislocations on different slip system planes. Single-crystal verification experiments are used to calibrate the CDD crystal plasticity model, and a set of simplified polycrystalline simulations demonstrates the model’s ability to capture the stress response from tensile experiments on α-iron.

## 1 Introduction

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.

## 2 Model Details

**F**

^{p}, while the elastic deformation gradient tensor,

**F**

^{e}, accounts for recoverable elastic stretch and rotation of the crystal lattice. The evolution of the plastic deformation is given as

**L**

^{p}is the plastic velocity gradient. The plastic velocity gradient is defined as the sum of the slip increments resulting from dislocation motion on all the slip systems [18]:

**s**

_{o}and

**m**

_{o}are the slip direction and slip plane normal unit vectors, respectively, in the reference configuration. The slip rate is connected to the glide of mobile dislocations via Orowan’s relation [19]:

*b*is the Burgers vector of the crystal, and $vglide\alpha $ is the glide velocity of the mobile dislocations. The evolution of the mobile dislocation rate is a constitutive expression described in Sec. 2.1. In this work, we apply a power law expression [20] for the dislocation glide velocity:

*v*

_{o}is the initial dislocation velocity,

*τ*

^{α}is the applied resolved shear stress on each slip system α, and

*g*

^{α}is the slip system resistance or strength, further defined in Eq. (6). The second Piola–Kirchoff stress tensor is used to calculate the applied shear stress, in a manner consistent with the use of the slip direction and slip plane normal unit vectors in Eq. (3).

*g*

^{α}in Eq. (5), provides a measure of the strength of each slip system to resist dislocation motion. Physically based frameworks write the constitutive slip system resistance equation as a function of dislocation and defect densities within the crystal [21,22]; both defects and other dislocations act as barriers to dislocation motion [23,24]. In this model, the resistance of the slip systems to dislocation motion is considered as the additive sum of the physical barriers to dislocation motion, including dislocation forests, and irradiation defects:

*g*

_{ps}is the isotropic Peierls strength of the material. For the lower temperature case in which only the {110} and {112} planes are active [26,27], we consider only these two slip plane types in our crystal plasticity model.

_{mbh}is a fitting parameter,

*μ*is the shear modulus of the material, and Ω is the interaction matrix containing the self- and latent-hardening parameters.

### 2.1 Dislocation Evolution.

*R*

_{c}, which is defined as a multiple of the Burgers vector.

*v*

_{glide}is the dislocation glide velocity from Eq. (5),

*R*

_{c}is the radius of capture for annihilating dislocations of opposite sign, and

*l*

_{inv}is the inverse mean free glide path. The inverse mean glide path is a measure of the distance between existing dislocations through which a mobile dislocation can glide [21]:

### 2.2 Dislocation Cross Slip Models.

*τ** is a critical stress for a dislocation to cross slip (in mm),

*τ*

^{(β)}is the applied shear stress on a slip system β (in mm),

*V*

_{a}is the volume (mm

^{3}) required for a dislocation to cross slip,

*k*is the Boltzmann constant, and

*T*is the temperature (in K). The probability of cross slip to move to a slip system is calculated for each slip system within a cross-slip family. Cross slip systems with higher applied stresses are more likely to receive cross slip dislocations.

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.”

#### 2.2.1 Stochastic-Only-Patra-McDowell Cross Slip Approach.

_{family}is the total number of slip systems within a cross slip family. The net change of cross slipped dislocations for a single slip system depends on two values: that of the probability of cross slip, and that of the mobile dislocation densities.

#### 2.2.2 Combined Stochastic-Monte Carlo Cross Slip Approach.

_{family}is the total number of slip systems within a cross slip family, $\rho received(\beta )$ is cross slip dislocations received from other slip systems in the same cross slip family, and $\rho given(\alpha )$ is the value of dislocation density that cross slips from the current (α) system to another system. $\rho received(\beta )$ is only nonzero if the Monte Carlo analysis determines that the current slip system (α) receives cross slip dislocations from a (β) system. If the random number generated its corresponding CDF bin are associated with the same slip system, then that slip system is considered to not produce any cross slipped dislocations, and $\rho given(\alpha )$ is zero.

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.

### 2.3 Irradiation Defect Evolution.

_{SIA}is a fitting parameter for SIA loop annihilation,

*R*

_{SIA}is the radius within which an SIA loop can capture and absorb dislocations, and

*ρ*

_{SIA}is the density of SIA loops. After an initial exposure to radiation, the density of SIA loops will decrease under mechanical loading of the crystal, as the dislocations absorb the SIA loops during the glide motion.

_{SIA}is a fitting coefficient. Here, we account for the impact of the 3D SIA loops by multiplying the density,

*ρ*

_{SIA}, by the average diameter,

*d*

_{SIA}, on the 2D glide slip planes of the mobile dislocations.

## 3 Model Calibration

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 (α_{5} in Eq. (15)) by varying the number of active slip systems: [001] has four, $[01\xaf1\xaf]$ has two, and $[3\xaf48]$ has only a single favored active slip system.

The single-crystal calibration simulations were performed on a 1 mm^{3} 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\xaf48]$ loading direction demonstrated the largest mesh dependence. The $[01\xaf1\xaf]$ and $[3\xaf48]$ 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^{−4} 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^{5} mm^{−2} [35] and distributed evenly across all slip systems.

Parameter | Value | Description |
---|---|---|

C_{11} | 242 × 10^{3} MPa | Elastic constant [23] |

C_{12} | 150 × 10^{3} MPa | Elastic constant [23] |

C_{44} | 112 × 10^{3} MPa | Elastic constant [23] |

μ | 80 × 10^{3} MPa | Shear modulus [23] |

b | 2.48 × 10^{−7} mm | Burgers vector [23] |

$\gamma \u02d9o$ | 4.0 × 10^{−2} | Reference strain rate |

m | 0.012 | Strain rate exponent |

α_{mbh} | 0.4 | Dispersed barrier |

Ω^{αα} | 1.0 | Self-hardening |

Ω^{αβ} | 0.2 | Latent-hardening |

R_{c} | 15b mm | Radius of capture |

τ* | 4 × 10^{−3} · μ MPa | Critical cross slip stress |

V_{a} | 20b^{3} mm^{3} | Cross slip volume |

k | 1.38065 × 10^{−20} MPa-mm^{3}/K | Bolztmann constant |

T | 298 K | Temperature |

Parameter | Value | Description |
---|---|---|

C_{11} | 242 × 10^{3} MPa | Elastic constant [23] |

C_{12} | 150 × 10^{3} MPa | Elastic constant [23] |

C_{44} | 112 × 10^{3} MPa | Elastic constant [23] |

μ | 80 × 10^{3} MPa | Shear modulus [23] |

b | 2.48 × 10^{−7} mm | Burgers vector [23] |

$\gamma \u02d9o$ | 4.0 × 10^{−2} | Reference strain rate |

m | 0.012 | Strain rate exponent |

α_{mbh} | 0.4 | Dispersed barrier |

Ω^{αα} | 1.0 | Self-hardening |

Ω^{αβ} | 0.2 | Latent-hardening |

R_{c} | 15b mm | Radius of capture |

τ* | 4 × 10^{−3} · μ MPa | Critical cross slip stress |

V_{a} | 20b^{3} mm^{3} | Cross slip volume |

k | 1.38065 × 10^{−20} MPa-mm^{3}/K | Bolztmann constant |

T | 298 K | Temperature |

The values of the dislocation mechanism terms, α_{1}–α_{6}, 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].

Loading direction | |||
---|---|---|---|

[100] | $[01\xaf1\xaf]$ | $[3\xaf48]$ | |

α_{1} | 0.03 | 0.03 | 0.03 |

α_{2} | 0.5 | 0.5 | 0.5 |

α_{3} | 0.002 | 0.002 | 0.002 |

α_{4} | 0.002 | 0.002 | 0.002 |

α_{5} | 0.015 | 0.0335 | 0.044 |

α_{6} | 1.0 | 1.0 | 1.0 |

g_{o} | 15 MPa | 8.8 MPa | 8.8 MPa |

Loading direction | |||
---|---|---|---|

[100] | $[01\xaf1\xaf]$ | $[3\xaf48]$ | |

α_{1} | 0.03 | 0.03 | 0.03 |

α_{2} | 0.5 | 0.5 | 0.5 |

α_{3} | 0.002 | 0.002 | 0.002 |

α_{4} | 0.002 | 0.002 | 0.002 |

α_{5} | 0.015 | 0.0335 | 0.044 |

α_{6} | 1.0 | 1.0 | 1.0 |

g_{o} | 15 MPa | 8.8 MPa | 8.8 MPa |

## 4 Discussion and Results

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.

### 4.1 Combined Stochastic-Monte Carlo and Stochastic-Only-Patra-McDowell Comparison.

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\xaf48]$ 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\xaf48]$ loading direction. In the stochastic-PM cross slip approach, we adjust the leading coefficient, α_{5}, 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\xaf11)[111]$ system—as well as, to a lesser extent, away from the secondary activated $(1\xaf21)[111\xaf]$ and $(1\xaf21)[111\xaf]$ 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\xaf11)[111]$ system mobile dislocation density. The saturation from the deterministic approach reduces the mobile dislocation density of the $(2\xaf11)[111]$ system to such an extent that the mobile dislocation density of the secondary $(1\xaf21)[111\xaf]$ system grows to parity with that of the primary $(2\xaf11)[111]$ system. This saturation of the $(2\xaf11)[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\xaf11)[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\xaf11)[111]$ system to grow at a moderate pace. Even though the mobile dislocation density on additional slip systems, $(1\xaf21)[111\xaf]$ and $(11\xaf2)[1\xaf11]$, increases, the favored slip system, $(2\xaf11)[111]$, consistently maintains the highest dislocation density value, as is expected in a single slip loading orientation.

### 4.2 Connection Between Cross Slip and Anisotropy.

Inclusion of anisotropy is necessary to capture the single slip system behavior of the $[3\xaf48]$ loading direction orientation. Similar cross slip probabilities are calculated on both the $(2\xaf11)[111]$ system and the $(1\xaf01)[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\xaf11)[111]$ to $(1\xaf01)[111]$ and $(1\xaf10)[111]$, as shown in Fig. 5(a). Transfer of dislocations from the active $(2\xaf11)[111]$ system to the $(1\xaf01)[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\xaf48]$ loading.

The higher slip system resistance allows the $(1\xaf01)[111]$ and $(1\xaf10)[111]$ systems to absorb the cross slip dislocations from the $(2\xaf11)[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\xaf48]$ direction.

### 4.3 Self-Interstitial Atom Loop Evolution Verification.

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^{3} 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^{13} mm^{−3}, 8.15 × 10^{12} mm^{−3}, 3.61 × 10^{12} mm^{−3}, 1.63 × 10^{12} mm^{−3}, 8.15 × 10^{11} mm^{−3}, and none (unirradiated case). As per Barton et al. [16], an initial dislocation density of 2 × 10^{−7} mm^{−2} 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.

### 4.4 Polycrystalline Simulations of Irradiated α-Iron.

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^{−6} 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^{7} mm^{−2}, which falls within the experimentally measured 7 ± 2 × 10^{7} mm^{−2} 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^{−4} 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.

## 5 Conclusions

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.

## Acknowledgment

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.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

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