## 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

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 [1417]. 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

As done by Asaro [18], our crystal plasticity model uses the multiplicative decomposition of the deformation gradient into elastic and plastic components:
$F=FeFp$
(1)
The crystal shape alteration due to the dislocation motion is accounted for by the plastic deformation gradient tensor, Fp, while the elastic deformation gradient tensor, Fe, accounts for recoverable elastic stretch and rotation of the crystal lattice. The evolution of the plastic deformation is given as
$F˙p=LpFp$
(2)
where Lp 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]:
$Lp=∑α=1nγ˙αsoα⊗moα$
(3)
where $γ˙α$ is the slip rate due to dislocation glide and so and mo 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]:
$γ˙α=ρmobileαbvglideα$
(4)
where $ρmobileα$ is the mobile dislocation density on each slip system α, b is the Burgers vector of the crystal, and $vglideα$ 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:
$vglideα=vob|ταgα|1/msgn(τα)ifτα≥goα$
(5)
where vo 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).
The slip system resistance, 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(α)=go(α)+gforest(α)+girradiation−defects(α)$
(6)
In Eq. (6), the initial lattice friction is represented by a constant, $goα$. The effect of anisotropy is accounted for by adjusting the initial slip system resistance for each type of slip system type [25]:
$go(α)={2.5⋅gpsgo{110}type1.0⋅gpsgo{112}type2.0⋅gpsgo{123}type$
(7)
where gps 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.
The slip system resistance due to dislocation-forest-type accumulation on all slip systems is modeled with a modified Bailey–Hirsch approach [21]:
$gforest(α)=αmbhbμ∑βΩ(αβ)(ρmobile(β)+ρimmobile(β))$
(8)
where the hardening is considered self-hardening when α = β and latent-hardening when α ≠ β. The coefficient α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.

Our crystal plasticity model is based on the CDD framework [28,29],with separate terms used to describe each of the specific physical interaction mechanisms in the dislocation evolution rate terms. The mobile dislocation evolution rate is fully coupled to the immobile dislocation evolution rate. The mobile dislocation evolution is governed by six terms. Each term in the equation represents a specific physical dislocation–interaction mechanism.
$ρ˙mobile(α)=ρ˙generation(α)−ρ˙mobile−annihilation(α)−ρ˙locking(α)+ρ˙freed(α)+ρ˙cross−slip(α)−ρ˙immobile−annihilation(α)$
(9)
The mobile dislocation evolution equation is coupled directly to the immobile dislocation evolution rate via three of the dislocation mechanism terms:
$ρ˙immobile(α)=ρ˙locking(α)−ρ˙freed(α)−ρ˙immobile−annihilation(α)$
(10)
The primary generation term for the mobile dislocation evolution rate, Eq. (9), is the creation of new mobile dislocations from Frank-Read sources:
$ρ˙generation(α)=α1ρmobile(α)vglidelinv$
(11)
The annihilation of two mobile dislocations of opposite sign is the second largest contributor to the mobile dislocation evolution rate. The annihilation interaction is allowed to occur only within a capture radius, Rc, which is defined as a multiple of the Burgers vector.
$ρ˙mobile−annihilation(α)=α2(ρmobile(α))22Rcbvglide$
(12)
The locking dislocations term, which is negative in the mobile dislocation evolution (Eq. (9)) acts as the source term for the immobile dislocation evolution, Eq. (10). This term represents the mobile dislocations that are then prevented from further movement by the dislocation locks.
$ρ˙locking(α)=α3ρmobile(α)vglidelinv$
(13)
The mobile and immobile dislocation evolution rates are also coupled via the freed dislocation mechanism term. This term represents those locked dislocations that are freed through the application of increasing shear stress on the slip system, and it transfers dislocations from the immobile to the mobile densities. The increasing shear stress supplies the additional driving force needed to overcome the strength of the locked barrier.
$ρ˙freed(α)=α4ρimm(α)|τ(α)g(α)|2vglidelinv$
(14)
The cross-slip term represents the movement of mobile dislocations from one slip system to another, within the same cross slip family. Depending on the applied shear stress in the receiving slip system, these cross slip dislocations may contribute to the mobile or immobile dislocation density.
$ρ˙cross−slip(α)=α5ρcross−slipped(α)vglidelinv$
(15)
The final mechanism term represents the annihilation between a mobile dislocation and an immobile dislocation of opposite sign. Unlike the mobile–mobile dislocation annihilation in Eq. (12), this term contributes to the coupling between the mobile and immobile dislocation evolution rates by removing equal quantities from both dislocation density populations.
$ρ˙immobile−annihilation(α)=α6ρmobile(α)ρimm(α)2Rcbvglide$
(16)
In these dislocation evolution expressions, the leading α coefficients are fitting parameters with values taken from dislocation dynamics simulation results [28], vglide is the dislocation glide velocity from Eq. (5), Rc is the radius of capture for annihilating dislocations of opposite sign, and linv 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]:
$linv(α)=βforest∑β(ρmobile(β)+ρimmobile(β))$
(17)
This mean free glide path quantity represents an inhibiting force on the dislocation motion in the evolution expressions in Eqs. (11), (13), (14), and (15), representing the increased difficulty of mobile dislocation glide as the total number of dislocations within a grain accumulate.

### 2.2 Dislocation Cross Slip Models.

Cross slip is a significant component of dislocation motion in BCC materials [30]. In lower length scale models such as dislocation dynamics, the probability of cross slip from single slip systems within a cross-slip family is calculated using a stochastic approach. Rhee et al. calculated the probability of cross slip as a function of applied shear stress and temperature [31].
$Pαβ=exp(−(τ*−|τ(β)|)VakT)$
(18)
where τ* is a critical stress for a dislocation to cross slip (in mm), τ(β) is the applied shear stress on a slip system β (in mm), Va is the volume (mm3) 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.

The stochastic-only cross slip approach is included in this work for comparison with the stochastic Monte Carlo combination approach. The stochastic-only model was implemented following the approach of Patra and McDowell [15]. In this approach, the probability of cross slip is also calculated using Eq. (18). The probability values are used directly as the coefficients for the mobile dislocation densities in order to calculate the cross slip dislocation density:
$ρcross−slipped(α)=Pαβ∑β≠αβ|familyρmobile(β)−ρmobile(α)∑β≠αβ|familyPαβ$
(19)
where α represents the current slip system, β represents the other slip systems in the same cross slip family, and β|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.

In the stochastic Monte Carlo combination approach to modeling dislocation cross slip, we calculate the probability of cross slip via the same function used by Rhee et al. Eq. (18) [31]. To perform the Monte-Carlo-type analysis, we construct a continuous distribution function (CDF) for each cross slip system by normalizing the probability by the sum of all probabilities within the slip system, such that the CDF is bounded between 0 and 1. A larger CDF bin indicates a correspondingly larger probability that a certain system will receive cross slipped dislocations. In the Monte Carlo analysis, performed for each slip system α in a cross slip family, a randomly generated number is compared against the bins in the CDF. The bin in which the random number falls is considered the slip system β that receives the cross slipped dislocations from slip system α:
$ρcross−slipped(α)=∑β=1β|familyPβαρ(β)=∑β=1β|familyρreceived(β)−ρgiven(α)$
(20)
where β|family is the total number of slip systems within a cross slip family, $ρreceived(β)$ is cross slip dislocations received from other slip systems in the same cross slip family, and $ρgiven(α)$ is the value of dislocation density that cross slips from the current (α) system to another system. $ρreceived(β)$ 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 $ρgiven(α)$ 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.

Irradiation defects contribute to the embrittlement of the slip systems by raising the slip system resistance, denoted as $girradiation−defects(α)$ in Eq. (6). In α-iron, SIA loops act as barriers to dislocation motion, and the 3D nature of these loops is accounted for in the barrier model with the cube root term:
$girradiation−defects(α)=αSIAμ(b2dSIAlSIA3)1/2$
(21)
where
$lSIA=1ρSIA[3]$
(22)
SIA loops can also be absorbed by dislocations that pass through the same plane on which an SIA loop exists. Following Mastorakos et al., the evolution of SIA loops is considered via a simple interaction model similar to the dislocation annihilation terms in Eqs. (12) and (16) [32]:
$ρ˙SIA=−4βSIARSIAρSIA∑αρmobile(α)vglide(α)$
(23)
where βSIA is a fitting parameter for SIA loop annihilation, RSIA 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.
As with forest dislocations, the SIA loops contribute to the shortening of the mean free glide path for mobile dislocations. We account for the contribution of the SIA loops to the mean free glide path by modifying Eq. (17) to be a function of both forest dislocations and SIA loops.
$linv(α)=βforest∑β(ρmobile(β)+ρimmobile(β))+βSIAρSIAdSIA$
(24)
As in Eq. (17), β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, dSIA, 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¯1¯]$ has two, and $[3¯48]$ has only a single favored active slip system.

The single-crystal calibration simulations were performed on a 1 mm3 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−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 × 105 mm−2 [35] and distributed evenly across all slip systems.

Fig. 1
Fig. 1
Close modal
Table 1

Values of the parameters used in the crystal plasticity model for single-crystal α-iron simulations

ParameterValueDescription
C11242 × 103 MPaElastic constant [23]
C12150 × 103 MPaElastic constant [23]
C44112 × 103 MPaElastic constant [23]
μ80 × 103 MPaShear modulus [23]
b2.48 × 10−7 mmBurgers vector [23]
$γ˙o$4.0 × 10−2Reference strain rate
m0.012Strain rate exponent
αmbh0.4Dispersed barrier
Ωαα1.0Self-hardening
Ωαβ0.2Latent-hardening
τ*4 × 10−3 · μ MPaCritical cross slip stress
Va20b3 mm3Cross slip volume
k1.38065 × 10−20 MPa-mm3/KBolztmann constant
T298 KTemperature
ParameterValueDescription
C11242 × 103 MPaElastic constant [23]
C12150 × 103 MPaElastic constant [23]
C44112 × 103 MPaElastic constant [23]
μ80 × 103 MPaShear modulus [23]
b2.48 × 10−7 mmBurgers vector [23]
$γ˙o$4.0 × 10−2Reference strain rate
m0.012Strain rate exponent
αmbh0.4Dispersed barrier
Ωαα1.0Self-hardening
Ωαβ0.2Latent-hardening
τ*4 × 10−3 · μ MPaCritical cross slip stress
Va20b3 mm3Cross slip volume
k1.38065 × 10−20 MPa-mm3/KBolztmann constant
T298 KTemperature

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

Table 2

Dislocation evolution parameters for Eqs. (9) and (10), calibrated for the crystal plasticity model in single-crystal α-iron simulations for three different loading directions

[100]$[01¯1¯]$$[3¯48]$
α10.030.030.03
α20.50.50.5
α30.0020.0020.002
α40.0020.0020.002
α50.0150.03350.044
α61.01.01.0
go15 MPa8.8 MPa8.8 MPa
[100]$[01¯1¯]$$[3¯48]$
α10.030.030.03
α20.50.50.5
α30.0020.0020.002
α40.0020.0020.002
α50.0150.03350.044
α61.01.01.0
go15 MPa8.8 MPa8.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¯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, α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¯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.

Fig. 2
Fig. 2
Close modal
Fig. 3
Fig. 3
Close modal

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.

Fig. 4
Fig. 4
Close modal

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.

### 4.2 Connection Between Cross Slip and Anisotropy.

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.

Fig. 5
Fig. 5
Close modal

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.

### 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 μm3 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 × 1013 mm−3, 8.15 × 1012 mm−3, 3.61 × 1012 mm−3, 1.63 × 1012 mm−3, 8.15 × 1011 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.

Table 3

Values of the parameters used in the CDD crystal plasticity model in the SIA loop terms, Eqs. (23) and (21), for verification of the SIA loop following Barton et al. [16]

ParameterValueDescription
αSIA0.7Hardening coeff. [38]
βSIA1.0Annihilation coeff. [32]
dSIA2.48 × 10−5 mmAverage diameter [16]
ParameterValueDescription
αSIA0.7Hardening coeff. [38]
βSIA1.0Annihilation coeff. [32]
dSIA2.48 × 10−5 mmAverage diameter [16]

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.

Fig. 6
Fig. 6
Close modal

### 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 × 107 mm−2, which falls within the experimentally measured 7 ± 2 × 107 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.

Table 4

Values of the parameters used in the SIA loop terms of the crystal plasticity model for polycrystalline α-iron exposed to a radiation dose of 0.1 dpa, including initial conditions

ParameterValueDescription
αSIA0.6Hardening coeff.
βSIA0.01Loop annihilation coeff.
ρSIA1.2 × 1012 mm−3Initial density [41]
dSIA7.0 × 10−6 mmAverage diameter [41]
ParameterValueDescription
αSIA0.6Hardening coeff.
βSIA0.01Loop annihilation coeff.
ρSIA1.2 × 1012 mm−3Initial density [41]
dSIA7.0 × 10−6 mmAverage diameter [41]

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

Fig. 7
Fig. 7
Close modal

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.

## References

1.
Allen
,
T. R.
, and
Busby
,
J. T.
,
2009
, “
Radiation Damage Concerns for Extended Light Water Reactor Service
,”
JOM
,
61
(
7
), pp.
29
34
.
2.
Zinkle
,
S. J.
, and
Ghoniem
,
N. M.
,
2011
, “
Prospects for Accelerated Development of High Performance Structural Materials
,”
J. Nucl. Mater.
,
417
(
1–3
), pp.
2
8
.
3.
Singh
,
B.
,
Foreman
,
A.
, and
Trinkaus
,
H.
,
1997
, “
,”
J. Nucl. Mater.
,
249
(
2–3
), pp.
103
115
.
4.
Zinkle
,
S. J.
, and
Matsukawa
,
Y.
,
2004
, “
Observation and Analysis of Defect Cluster Production and Interactions With Dislocations
,”
J. Nucl. Mater.
,
329–333
(
Part A
), pp.
88
96
.
5.
de la Rubia
,
T.
,
Zbib
,
H.
,
Khraishi
,
T.
,
Wirth
,
B.
,
Victoria
,
M.
, and
Caturla
,
M.
,
2000
, “
Multiscale Modelling of Plastic Flow Localization in Irradiated Materials
,”
Nature
,
406
(
6798
), pp.
871
874
.
6.
Klueh
,
R.
, and
Nelson
,
A.
,
2007
, “
Ferritic/Martensitic Steels for Next-Generation Reactors
,”
J. Nucl. Mater.
,
371
(
1
), pp.
37
52
. Nuclear Fuels and Structural Materials 1.
7.
Arsenlis
,
A.
, and
Parks
,
D. M.
,
2002
, “
Modeling the Evolution of Crystallographic Dislocation Density in Crystal Plasticity
,”
J. Mech. Phys. Solids
,
50
(
9
), pp.
1979
2009
.
8.
Cheong
,
K.-S.
, and
Busso
,
E. P.
,
2004
, “
Discrete Dislocation Density Modelling of Single Phase FCC Polycrystal Aggregates
,”
Acta Mater.
,
52
(
19
), pp.
5665
5675
.
9.
Kocks
,
U.
,
1976
, “
Laws for Work-Hardening and Low-Temperature Creep
,”
J. Eng. Mater.
,
98
(
1
), pp.
76
85
.
10.
Roters
,
F.
,
Raabe
,
D.
, and
Gottstein
,
G.
,
2000
, “
Work Hardening in Heterogeneous Alloys – a Microstructural Approach Based on Three Internal State Variables
,”
Acta Mater.
,
48
(
17
), pp.
4181
4189
.
11.
Ma
,
A.
, and
Roters
,
F.
,
2004
, “
A Constitutive Model for FCC Single Crystals Based on Dislocation Densities and Its Application to Uniaxial Compression of Aluminium Single Crystals
,”
Acta Mater.
,
52
(
12
), pp.
3603
3612
.
12.
Arsenlis
,
A.
,
Wirth
,
B. D.
, and
Rhee
,
M.
,
2004
, “
Dislocation Density-Based Constitutive Model for the Mechanical Behaviour of Irradiated Cu
,”
Philos. Mag.
,
84
(
34
), pp.
3617
3635
.
13.
Wirth
,
B. D.
,
Caturla
,
M. J.
,
Diaz de la Rubia
,
T.
,
Khraishi
,
T.
, and
Zbib
,
H.
,
2001
, “
,”
Nucl. Instrum. Methods Phys. Res. Sect. B: Beam Interact. Mater. Atoms
,
180
(
1–4
), pp.
23
31
.
14.
Vincent
,
L.
,
Libert
,
M.
,
Marini
,
B.
, and
Rey
,
C.
,
2010
, “
Towards a Modelling of RPV Steel Brittle Fracture Using Crystal Plasticity Computation on Polycrystalline Aggregates
,”
J. Nucl. Mater.
,
406
(
1
), pp.
91
96
.
15.
Patra
,
A.
, and
McDowell
,
D. L.
,
2013
, “
Continuum Modeling of Localized Deformation in Irradiated BCC Materials
,”
J. Nucl. Mater.
,
432
(
1
), pp.
414
427
.
16.
Barton
,
N. R.
,
Arsenlis
,
A.
, and
Marian
,
J.
,
2013
, “
A Polycrystal Plasticity Model of Strain Localization in Irradiated Iron
,”
J. Mech. Phys. Solids
,
61
(
2
), pp.
341
351
.
17.
Chakraborty
,
P.
, and
Biner
,
S. B.
,
2016
, “
Crystal Plasticity Modeling of Irradiation Effects on Flow Stress in Pure-Iron and Iron-Copper Alloys
,”
Mech. Mater.
,
101
, pp.
71
80
.
18.
Asaro
,
R. J.
,
1983
, “
Crystal Plasticity
,”
J. Appl. Mech.
,
50
(
4b
), pp.
921
934
.
19.
Orowan
,
E.
,
1940
, “
Problems of Plastic Gliding
,”
Proc. Phys. Soc.
,
52
(
1
), p.
8
.
20.
Asaro
,
R. J.
, and
Needleman
,
A.
,
1985
, “
Overview No. 42: Texture Development and Strain Hardening in Rate Dependent Polycrystals
,”
Acta Metall.
,
33
(
6
), pp.
923
953
.
21.
Ohashi
,
T.
,
1994
, “
Numerical Modelling of Plastic Multislip in Metal Crystals of FCC Type
,”
Philos. Mag. A
,
70
(
5
), pp.
793
803
.
22.
Roters
,
F.
,
Eisenlohr
,
P.
,
Hantcherli
,
L.
,
Tjahjanto
,
D. D.
,
Bieler
,
T. R.
, and
Raabe
,
D.
,
2010
, “
Overview of Constitutive Laws, Kinematics, Homogenization and Multiscale Methods in Crystal Plasticity Finite-Element Modeling: Theory, Experiments, Applications
,”
Acta Mater.
,
58
(
4
), pp.
1152
1211
.
23.
Hirth
,
J. P.
, and
Lothe
,
J.
,
1982
,
Theory of Dislocations
, 2nd ed.,
John Wiley & Sons
,
Malabar, FL
.
24.
Zbib
,
H. M.
,
Dıaz De La Rubia
,
T.
, and
Rhee
,
M.
,
2000
, “
3D Dislocation Dynamics: Stress-Strain Behavior and Hardening Mechanisms in FCC and BCC Metals
,”
J. Nucl. Mater.
,
276
(
1
), pp.
154
165
.
25.
Lyu
,
H.
,
Ruimi
,
A.
, and
Zbib
,
H. M.
,
2015
, “
A Dislocation-Based Model for Deformation and Size Effect in Multi-Phase Steels
,”
Int. J. Plast.
,
72
(
4
), pp.
44
59
.
26.
Raabe
,
D.
,
1995
, “
Simulation of Rolling Textures of BCC Metals Considering Grain Interactions and Crystallographic Slip on {110}, {112} and {123} Planes
,”
Mater. Sci. Eng. A
,
197
(
1
), pp.
31
37
.
27.
Alankar
,
A.
,
Field
,
D. P.
, and
Raabe
,
D.
,
2014
, “
Plastic Anisotropy of Electro-Deposited Pure α-Iron With Sharp Crystallographic < 1 1 1 > Texture in Normal Direction: Analysis by An Explicitly Dislocation-Based Crystal Plasticity Model
,”
Int. J. Plast.
,
52
(
Special Issue in Honor of Hussein Zbib
), pp.
18
32
.
28.
Li
,
D.
,
Zbib
,
H.
,
Sun
,
X.
, and
Khaleel
,
M.
,
2014
, “
Predicting Plastic Flow and Irradiation Hardening of Iron Single Crystal with Mechanism-Based Continuum Dislocation Dynamics
,”
Int. J. Plast.
,
52
(
8
), pp.
3
17
.
29.
Taheri-Nassaj
,
N.
, and
Zbib
,
H. M.
,
2016
, “
A Mesoscale Model of Plasticity: Dislocation Dynamics and Patterning (One-Dimensional)
,”
J. Eng. Mater.
,
138
(
4
), p.
041015
.
30.
Püschl
,
W.
,
2002
, “
Models for Dislocation Cross-Slip in Close-Packed Crystal Structures: a Critical Review
,”
Prog. Mater. Sci.
,
47
(
4
), pp.
415
461
.
31.
Rhee
,
M.
,
Zbib
,
H. M.
,
Hirth
,
J. P.
, and
Huang
,
H.
,
1998
, “
Models for Long-/Short-Range Interactions and Cross Slip in 3D Dislocation Simulation of BCC Single Crystals
,”
Modell. Simul. Mater. Sci. Eng.
,
6
(
4
), p.
467
.
32.
Mastorakos
,
I.
, and
Zbib
,
H.
,
2014
, “
A Multiscale Approach to Study the Effect of Chromium and Nickel Concentration in the Hardening of Iron Alloys
,”
J. Nucl. Mater.
,
449
(
1
), pp.
101
110
.
33.
Permann
,
C. J.
,
Gaston
,
D. R.
,
Andrš
,
D.
,
Carlsen
,
R. W.
,
Kong
,
F.
,
Lindsay
,
A. D.
,
Miller
,
J. M.
,
Peterson
,
J. W.
,
Slaughter
,
A. E.
,
Stogner
,
R. H.
, and
Martineau
,
R. C.
,
2020
, “
MOOSE: Enabling Massively Parallel Multiphysics Simulation
,”
SoftwareX
,
11
(
10
), p.
100430
.
34.
Keh
,
A.
,
1965
, “
Work Hardening and Deformation Sub-structure in Iron Single Crystals Deformed in Tension At 298 K
,”
Philos. Mag.
,
12
(
115
), pp.
9
30
.
35.
Lee
,
M.
,
Lim
,
H.
,
,
B.
,
Hirth
,
J.
, and
Wagoner
,
R.
,
2010
, “
A Dislocation Density-Based Single Crystal Constitutive Equation
,”
Int. J. Plast.
,
26
(
7
), pp.
925
938
.
36.
Pitts
,
S. A.
,
2019
,
Modeling and Simulation of Microstructure Evolution and Deformation in an Irradiated Environment, Publication Number 13810450, Doctoral dissertation
,
ProQuest Dissertations Publishing
,
Washington State University, Pullman, WA
.
37.
Arsenlis
,
A.
,
Rhee
,
M.
,
Hommes
,
G.
,
Cook
,
R.
, and
Marian
,
J.
,
2012
, “
A Dislocation Dynamics Study of the Transition From Homogeneous to Heterogeneous Deformation in Irradiated Body-Centered Cubic Iron
,”
Acta Mater.
,
60
(
9
), pp.
3748
3757
.
38.
Lambrecht
,
M.
,
Meslin
,
E.
,
Malerba
,
L.
,
Hernández-Mayoral
,
M.
,
Bergner
,
F.
,
Pareige
,
P.
,
,
B.
, and
Almazouzi
,
A.
,
2010
, “
On the Correlation Between Irradiation-Induced Microstructural Features and the Hardening of Reactor Pressure Vessel Steels
,”
J. Nucl. Mater.
,
406
(
1
), pp.
84
89
.
39.
Lambrecht
,
M.
,
Malerba
,
L.
, and
Almazouzi
,
A.
,
2008
, “
Influence of Different Chemical Elements on Irradiation-Induced Hardening Embrittlement of RPV Steels
,”
J. Nucl. Mater.
,
378
(
3
), pp.
282
290
.
40.
Python Software Foundation
,
2018
, “
The Python Standard Library
,” https://www.python.org/, Accessed July 18, 2018.
41.
Meslin
,
E.
,
Lambrecht
,
M.
,
Hernández-Mayoral
,
M.
,
Bergner
,
F.
,
Malerba
,
L.
,
Pareige
,
P.
,
,
B.
,
Barbu
,
A.
,
Gómez-Briceño
,
D.
,
Ulbricht
,
A.
, and
Almazouzi
,
A.
,
2010
, “
Characterization of Neutron-Irradiated Ferritic Model Alloys and a RPV Steel From Combined APT, SANS, TEM and PAS Analyses
,”
J. Nucl. Mater.
,
406
(
1
), pp.
73
83
.