Microstructural comparison of the kinematics of discrete and continuum dislocations models

From a geometrical viewpoint, a dislocation line is an oriented curve $\boldsymbol{c}(\ell )$ which can be parameterized by its arc-length $\ell$. The local unit tangent vectors are the first derivative of $\boldsymbol{c}$

$$\begin{eqnarray}&&\boldsymbol{l}^{{\text{d}}}(\ell )=\frac{\text{d}\boldsymbol{c}}{\text{d}\ell}\,,\end{eqnarray} \tag{ 1 }$$

while the curvature vector is given by the change of orientation per unit arc-length,

$$\begin{eqnarray}&&\boldsymbol{k}^{{\text{d}}}=\frac{\text{d}\boldsymbol{l}^{{\text{d}}}}{\text{d}\ell}=\frac{{{\text{d}}^{2}}\boldsymbol{c}}{\text{d}{{\ell}^{2}}}\,.\end{eqnarray} \tag{ 2 }$$

The scalar curvature is the reciprocal of the local radius of curvature R of the dislocation line:

$$\begin{eqnarray}&&{{k}^{\text{d}}}=\,\parallel \boldsymbol{k}^{{\text{d}}}\parallel \,=\frac{1}{R}.\end{eqnarray} \tag{ 3 }$$

Although parameterization by arc-length is the most natural choice for defining geometrical properties of a curve, in general the parameter $\ell$ is only available after the curve has been constructed. Thus, for practical purposes, it is often necessary to prescribe a curve through a parameter u, say $u\in [0,1]$. In this case unit tangent and curvature become

$$\begin{eqnarray}&&\boldsymbol{l}^{{\text{d}}}(u)=\frac{1}{J}\frac{\text{d}\boldsymbol{c}}{\text{d}u}\quad \text{and}\quad \boldsymbol{k}^{{\text{d}}}(u)=\frac{1}{{{J}^{2}}}\frac{{{\text{d}}^{2}}\boldsymbol{c}}{\text{d}{{u}^{2}}}-\frac{1}{{{J}^{3}}}\frac{\text{d}J}{\text{d}u}\frac{\text{d}\boldsymbol{c}}{\text{d}u}\end{eqnarray} \tag{ 4 }$$

where the scalar J and its derivative are given by

$$\begin{eqnarray}&&J(u)=\frac{\text{d}\ell}{\text{d}u}=\parallel \frac{\text{d}\boldsymbol{c}}{\text{d}u}\parallel \quad \text{and}\quad \frac{\text{d}J}{\text{d}u}=\frac{1}{J}\frac{\text{d}\boldsymbol{c}}{\text{d}u}\cdot \frac{{{\text{d}}^{2}}\boldsymbol{c}}{\text{d}{{u}^{2}}}\,.\end{eqnarray} \tag{ 5 }$$

In this work dislocation lines are discretized into segments, each of which is represented by a cubic Hermite spline of the form

$$\begin{eqnarray}\boldsymbol{c}(u)=\left[\begin{array}{c} 1 & u & {{u}^{2}} & {{u}^{3}} \\ \end{array}\right]\left[\begin{array}{c} 1 & 0 & 0 & 0 \\ 0 & \gamma & 0 & 0 \\ -3 & -2\gamma & 3 & -\gamma \\ 2 & \gamma & -2 & \gamma \\ \end{array}\right]\left[\begin{array}{*{35}{l}} \boldsymbol{P}_{{0}} \\ \boldsymbol{T}_{{0}} \\ \boldsymbol{P}_{{1}} \\ \boldsymbol{T}_{{1}} \end{array}\right],\end{eqnarray} \tag{ 6 }$$

where $\left(\boldsymbol{P}_{{0}},\boldsymbol{T}_{{0}}\right)$ are position and parametric tangent vectors of the first knot of the segment (u  =  0), $\left(\boldsymbol{P}_{{1}},\boldsymbol{T}_{{1}}\right)$ are the corresponding quantities for the second knot of the segment (u  =  1), $\gamma =\parallel \boldsymbol{P}_{{1}}-\boldsymbol{P}_{{0}}{{\parallel}^{\alpha}}$, and α is a tension parameter³. Cubic splines are used instead of linear segments in order to have non-vanishing curvature along each segment. This type of representation is adopted in the so-called Parametric Dislocation Dynamics method (PDD) [11]. The numerical studies below are based on the existing implementation [14, 46].

As an auxiliary functional for averaging geometrical properties of discrete dislocations, we introduce the Dirac delta-operator on a dislocation line $\boldsymbol{c}$ as

$$\begin{eqnarray}&&{{\delta _c}\left( {[{\rm{ \bullet }}],r} \right) = \mathop \int \nolimits _0^{{L_c}}\delta \left( {c(\ell ) - r} \right){\rm{ \bullet d}}\ell ,}\end{eqnarray} \tag{ 7 }$$

where the symbol '•' is a placeholder for an arbitrary line property described by a scalar, vectorial or tensorial function of $\boldsymbol{c}$ and r is a point in space. Using (7), for each discrete line $\boldsymbol{c}$ we can define the scalar dislocation density ${{\rho}_{\boldsymbol{c}}}$, the line direction density (or vector of signed GND densities) $\boldsymbol{\varrho }_{{\boldsymbol{c}}}$, the Kröner–Nye tensor $\boldsymbol{\alpha }_{{\boldsymbol{c}}}$ [47], and the curvature density vector $\boldsymbol{q}_{{\boldsymbol{c}}}$ as⁴

$$\begin{eqnarray}&&{{\rho}_{\boldsymbol{c}}}\left(\boldsymbol{r}\right)={{\delta}_{\boldsymbol{c}}}\left([1],\boldsymbol{r}\right)=\int_{0}^{{{L}_{\boldsymbol{c}}}}\delta \left(\boldsymbol{c}(\ell )-\boldsymbol{r}\right)\text{d}\ell \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\boldsymbol{\varrho }_{{\boldsymbol{c}}}\left(\boldsymbol{r}\right)={{\delta}_{\boldsymbol{c}}}\left(\left[\frac{\text{d}\boldsymbol{c}}{\text{d}\ell}\right],\boldsymbol{r}\right)=\int_{0}^{{{L}_{\boldsymbol{c}}}}\delta \left(\boldsymbol{c}(\ell )-\boldsymbol{r}\right)\frac{\text{d}\boldsymbol{c}}{\text{d}\ell}\,\text{d}\ell \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\boldsymbol{\alpha }_{{\boldsymbol{c}}}\left(\boldsymbol{r}\right)={{\delta}_{\boldsymbol{c}}}\left(\left[\frac{\text{d}\boldsymbol{c}}{\text{d}\ell}\otimes \boldsymbol{b}\right],\boldsymbol{r}\right)=\int_{0}^{{{L}_{\boldsymbol{c}}}}\delta \left(\boldsymbol{c}(\ell )-\boldsymbol{r}\right)\frac{\text{d}\boldsymbol{c}}{\text{d}\ell}\otimes \boldsymbol{b}\,\text{d}\ell \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&\boldsymbol{q}_{{\boldsymbol{c}}}\left(\boldsymbol{r}\right)={{\delta}_{\boldsymbol{c}}}\left(\left[\frac{{{\text{d}}^{2}}\boldsymbol{c}}{\text{d}{{\ell}^{2}}}\right],\boldsymbol{r}\right)=\int_{0}^{{{L}_{\boldsymbol{c}}}}\delta \left(\boldsymbol{c}(\ell )-\boldsymbol{r}\right)\,\frac{{{\text{d}}^{2}}\boldsymbol{c}}{\text{d}{{\ell}^{2}}}\,\text{d}\ell .\end{eqnarray} \tag{ 11 }$$

The later on used scalar, signed curvature density ${{q}_{\boldsymbol{c}}}$ can be obtained by projecting the curvature vector on the outwards pointing normal

$$\begin{eqnarray}&&{{q}_{\boldsymbol{c}}}\left(\boldsymbol{r}\right)=\int_{0}^{{{L}_{\boldsymbol{c}}}}\delta \left(\boldsymbol{c}(\ell )-\boldsymbol{r}\right)\,\frac{{{\text{d}}^{2}}\boldsymbol{c}}{\text{d}{{\ell}^{2}}}\cdot \left(\frac{\text{d}\boldsymbol{c}}{\text{d}\ell}\times \boldsymbol{n}\right)\,\text{d}\ell ,\end{eqnarray} \tag{ 12 }$$

where $\boldsymbol{n}$ is the slip plane normal. These discrete measures are well suited for spatial (or ensemble) averaging over a number of lines, and we introduce an averaging operator for the volume ${{V}_{\boldsymbol{r}}}$ of size V as

$$\begin{eqnarray}&&{{\langle \bullet \rangle}_{V,\boldsymbol{r}}}:=(1/V){{\mathop{\int}^{}}_{{{V}_{\boldsymbol{r}}}}}\bullet\;{{\text{d}}^{3}}r,\end{eqnarray} \tag{ 13 }$$

where the averaging volume is centered around $\boldsymbol{r}$. Application of (13) to (8)–(12) leads to the definition of the corresponding average fields:

$$\begin{eqnarray}&&\rho \left(\boldsymbol{r}\right)={{\langle \underset{\boldsymbol{c}}{\mathop \sum}\,{{\rho}_{\boldsymbol{c}}}\rangle}_{V,\boldsymbol{r}}}=\frac{1}{V}\underset{\boldsymbol{c}}{\mathop \sum}\,\int_{\mathcal{L}_{\boldsymbol{c}}^{V,\boldsymbol{r}}}\,\text{d}\ell \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\boldsymbol{\varrho }\left(\boldsymbol{r}\right)=\langle \underset{\boldsymbol{c}}{\mathop \sum}\,\boldsymbol{\varrho }_{{\boldsymbol{c}}}\rangle =\frac{1}{V}\underset{\boldsymbol{c}}{\mathop \sum}\,\int_{\mathcal{L}_{\boldsymbol{c}}^{V,\boldsymbol{r}}}\frac{\text{d}\boldsymbol{c}}{\text{d}\ell}\,\text{d}\ell \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\boldsymbol{\alpha }\left(\boldsymbol{r}\right)=\langle \underset{\boldsymbol{c}}{\mathop \sum}\,\boldsymbol{\alpha }_{{\boldsymbol{c}}}\rangle =\frac{1}{V}\underset{\boldsymbol{c}}{\mathop \sum}\,\int_{\mathcal{L}_{\boldsymbol{c}}^{V,\boldsymbol{r}}}\frac{\text{d}\boldsymbol{c}}{\text{d}\ell}\otimes \boldsymbol{b}\,\text{d}\ell \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&q\left(\boldsymbol{r}\right)=\langle \underset{\boldsymbol{c}}{\mathop \sum}\,{{q}_{\boldsymbol{c}}}\rangle =\frac{1}{V}\underset{\boldsymbol{c}}{\mathop \sum}\,\int_{\mathcal{L}_{\boldsymbol{c}}^{V,\boldsymbol{r}}}\,\frac{{{\text{d}}^{2}}\boldsymbol{c}}{\text{d}{{\ell}^{2}}}\cdot \left(\frac{\text{d}\boldsymbol{c}}{\text{d}\ell}\times \boldsymbol{n}\right)\,\text{d}\ell \end{eqnarray} \tag{ 17 }$$

Therein, $\mathcal{L}_{\boldsymbol{c}}^{V,\boldsymbol{r}}\subset \boldsymbol{c}$ denotes a section of line $\boldsymbol{c}$ contained inside the averaging volume V which is centered at $\boldsymbol{r}$. Note, that upon volume integration of the total density (14) the total dislocation line length contained within the volume is obtained. Analogously, integration of the signed curvature density (17) yields the number of closed loops as multiple of $2\pi$. The vector of signed GND densities $\boldsymbol{\varrho }$ contains densities of positive and negative edge and screw dislocations as later used in the CDD evolution equation (23). We remark that the (scalar) geometrically necessary dislocation density ${{\rho}_{\text{G}}}$ is the norm of $\boldsymbol{\varrho }$, i.e. ${{\rho}_{\text{G}}}\equiv \varrho =\,\parallel \boldsymbol{\varrho }\parallel$, and that the average line direction of the geometrically necessary dislocations is given by the unit vector $\boldsymbol{l}_{{\boldsymbol{\varrho }}}=\boldsymbol{\varrho }/\parallel \boldsymbol{\varrho }\parallel$.

Numerically, the fields (14)–(17) can be discretized by replacing the spatial points $\boldsymbol{r}$ with a finite set of points $\boldsymbol{r}_{{i}}$ which have the coordinates $\left({{x}_{i}},{{y}_{i}},{{z}_{i}}\right)$. On a regular grid each point $\boldsymbol{r}_{{i}}$ is then the center of a sub-volume ${{V}_{i}}=\Delta x \Delta y \Delta z$ which defines the averaging sub-domain ${{\Omega }_{i}}$:

$$\begin{eqnarray}&&{{\Omega }_{i}}=\left[{{x}_{i}}-\frac{1}{2}\Delta x,{{x}_{i}}+\frac{1}{2}\Delta x\right]\times \ldots \times \left[{{z}_{i}}-\frac{1}{2}\Delta z,{{z}_{i}}+\frac{1}{2}\Delta x\right]\end{eqnarray} \tag{ 18 }$$

and determines the resolution of the continuum representation. As a consequence of this coarse graining, it is well possible that inside a sub-domain ${{\Omega }_{i}}$ line segments 'cross each other' without intersections: although these lines are located on different physical glide planes they can not be resolved separately after coarse graining⁵. The discretized version of the scalar dislocation density (14) is:

$$\begin{eqnarray}&&{{\rho}_{i}}=\rho \left(\boldsymbol{r}_{{i}}\right)=\frac{1}{{{V}_{i}}}\underset{\boldsymbol{c}}{\mathop \sum}\,\left(\underset{k,\,\boldsymbol{c}\left({{u}_{k}}\right)\in {{\Omega }_{i}}}{\mathop \sum}\,{{J}_{\boldsymbol{c}}}\left({{u}_{k}}\right){{w}_{k}}\right).\end{eqnarray} \tag{ 19 }$$

In (19) we compute the line length contribution of each dislocation segment by numerical quadrature: letting u_k and w_k be the abscissas and weights of the quadrature method of choice⁶ in the unit interval [0, 1], respectively, each quadrature point contributes a line length $J\left({{u}_{k}}\right){{w}_{k}}$ to the domain ${{\Omega }_{i}}$ which contains $\boldsymbol{c}\left({{u}_{k}}\right)$. Conveniently, the quantity J(u_k) can be directly computed for each dislocation segment from (5) and (6). Analogously, we define the discretized vector of signed GND densities, and the discretized signed curvature density as

$$\begin{eqnarray}&&\boldsymbol{\varrho }_{{i}}=\boldsymbol{\varrho }\left(\boldsymbol{r}_{{i}}\right)=\frac{1}{{{V}_{i}}}\underset{\boldsymbol{c}}{\mathop \sum}\,\left(\underset{k,\,\boldsymbol{c}\left({{u}_{k}}\right)\in {{\Omega }_{i}}}{\mathop \sum}\,\frac{\text{d}\boldsymbol{c}}{\text{d}\ell}\left({{u}_{k}}\right)\,{{J}_{\boldsymbol{c}}}\left({{u}_{k}}\right){{w}_{k}}\right)\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{{q}_{i}}=q\left(\boldsymbol{r}_{{i}}\right)=\frac{1}{{{V}_{i}}}\underset{\boldsymbol{c}}{\mathop \sum}\,\left(\underset{k,\,\boldsymbol{c}\left({{u}_{k}}\right)\in {{\Omega }_{i}}}{\mathop \sum}\,{{\left[\frac{{{\text{d}}^{2}}\boldsymbol{c}}{\text{d}{{\ell}^{2}}}\cdot \left(\frac{\text{d}\boldsymbol{c}}{\text{d}\ell}\times \boldsymbol{n}\right)\right]}_{{{u}_{k}}}}{{J}_{\boldsymbol{c}}}\left({{u}_{k}}\right){{w}_{k}}\right).\end{eqnarray} \tag{ 21 }$$

The Continuum Dislocation Dynamics theory provides conceptual steps to map a discrete dislocation system onto a set of continuous, density-like field variables together with the respective evolution equations governing their change in time. The original, higher-dimensional theory (hdCDD) was defined in a configuration space which added an additional dimension to the spatial domain [32, 48]. hdCDD is very accurate and contains very detailed microstructure information [39, 49] but suffers from requiring a large number of computational degrees of freedoms. As a computationally favorable albeit somewhat less accurate theory a simplified version (CDD) is used, which is based on spatial evolution equations for the total dislocation density $\rho$, a vector of geometrically necessary signed densities $\boldsymbol{\varrho }$ and the curvature density $q$ [35, 45, 50]. We assume that the components of $\boldsymbol{\varrho }$ are oriented such that they represent the orientation of screw and edge dislocations: $\boldsymbol{\varrho }=\left[{{\varrho}_{\text{s}}},{{\varrho}_{\text{e}}}\right]$. The temporal evolution of these field quantities are—in local slip system coordinates—governed by the following set of equations:

$$\begin{eqnarray}&&{{\partial}_{t}}\rho =-\nabla \cdot \left(v\boldsymbol{\varrho }^{{\bot}}\right)+vq\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{{\partial}_{t}}\boldsymbol{\varrho }=-\nabla \times \left(v\rho \boldsymbol{n}\right)\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{{\partial}_{t}}q=-\nabla \cdot \left(-v{{\mathbf{Q}}^{(1)}}+{{\mathbf{A}}^{(2)}}\cdot \nabla v\right),\end{eqnarray} \tag{ 24 }$$

where $\boldsymbol{n}$ denotes the slip plane normal. The by 90° rotated GND density vector is $\boldsymbol{\varrho }^{{\bot}}=\left[{{\varrho}_{\text{e}}},-{{\varrho}_{\text{s}}}\right]$, and we assume ${{\mathbf{Q}}^{(1)}}=-\boldsymbol{\varrho }^{{\bot}}q/\rho$ (see [50] for a discussion of this assumption). The evolution equations are mathematically closed by an expression for the tensor ${{\mathbf{A}}^{(2)}}$ which here is obtained from a 'maximum entropy' approach [50]:

$$\begin{eqnarray}&&{{\mathbf{A}}^{(2)}}=\frac{\rho}{2}\left[(1+\Phi)\boldsymbol{l}_{{\varrho}}\otimes \boldsymbol{l}_{{\varrho}}+(1-\Phi)\boldsymbol{l}_{{{{\varrho}^{\bot}}}}\otimes \boldsymbol{l}_{{{{\varrho}^{\bot}}}}\right]\end{eqnarray} \tag{ 25 }$$

where $\boldsymbol{l}_{{{{\varrho}^{\bot}}}}$ is the unit vector perpendicular to $\boldsymbol{l}_{{\varrho}}$. As shown by Monavari et al [50] $\Phi$ can be approximated as $\Phi\approx {{(\varrho /\rho )}^{2}}\left(1+{{(\varrho /\rho )}^{4}}\right)/2$ and is a non-linear interpolation between $\varrho =0$ (isotropic dislocation arrangement) and $\varrho =\rho$ (fully polarized dislocation arrangement). In section 6 we numerically investigate the validity of this approximation. The corresponding non-dimensionless scaling of the CDD equations (22)–(24) can be found in appendix D of [51].

In general, the velocity v in the above equations needs to be specified in terms of a constitutive equation through which elastic interactions enter the model. This is a priori not part of CDD: e.g. long-range (or 'internal') interaction stresses have to be obtained from the additional solution of a dislocation eigenstrain problem [42, 52, 53] while elastic short-range interactions have to be recovered in an alternative way based on CDD fields (see e.g. [54–56] for steps into this direction) due to the limited resolution of any continuous description. Within the present work we will not tackle this problem of 'dynamic closure' and assume for benchmark purposes velocity fields given as stationary analytical functions that do not depend on the dislocation state⁷. In other words: we ignore all dislocation interactions and concentrate exclusively on the 'kinematics', that is, how curved and connected lines move and evolve in space and time without interactions.

Equations (19)–(21) can be used to extract coarse grained fields of geometrical data from a given configuration of discrete dislocations. However, the CDD evolution equations contain spatial derivatives which have to be approximated numerically and which thus require field data with a certain level of spatial smoothness. It is therefore helpful to replace the Dirac delta-function $\delta \left(\boldsymbol{r}\right)$ in the convolution integral (7) with a regularization function $G\left(\boldsymbol{r}\right)$. To improve the numerical efficiency during computation of the convolution it is beneficial to make use of the sifting property of the Dirac function: we can always write the regularized function as $G=G*\delta$, where the symbol $*$ indicates convolution over three-dimensional space. Together with the fact that convolution and volume averaging commute, it follows that regularized fields can be obtained by convolution of (19)–(21) with G. In our numerical implementation, we chose a Gaussian standard distribution as the regularization function:

$$\begin{eqnarray}&&G\left(\boldsymbol{r}\right)=\frac{1}{s\sqrt{2\pi}}\exp \left(-\frac{\parallel \boldsymbol{r}{{\parallel}^{2}}}{2{{s}^{2}}}\right),\end{eqnarray} \tag{ 26 }$$

where the standard deviation s characterizes the width of the dislocation density distribution. For numerical reasons G is approximated by a discrete Gauss function ${{G}^{\text{d}}}$ which is non-zero only at a finite number of discrete points r_i and is normalized such that $\Delta x \Delta y \Delta z\sum_{i}{{G}^{\text{d}}}(i)={{\mathop{\int}^{}}_{x,y}}G\left(r\left({{x}^{\prime}},{{y}^{\prime}},{{z}^{\prime}}\right)\right)\text{d}{{x}^{\prime}}\,\text{d}{{y}^{\prime}}\,\text{d}{{z}^{\prime}}$ within a small numerical tolerance. We then compute the convolution of the coarse grained fields (equations (19)–(21)) with the discrete Gauss function ${{G}^{\text{d}}}$. We note, that in (26) the only free parameter is the standard deviation s which subsequently will be chosen as the average dislocation spacing $\bar{s}$ (the mean separation between a dislocation and its nearest neighbors) inside the averaging volume (see also the discussion in [53]). Note that choosing a much larger value for s would destroy details of the dislocation structure to a large extend while if s were very small this would result in a quasi-discrete rather than continuous density structure. This is further discussed and demonstrated in the appendix A.

The following examples study a system as sketched in figure 1 where a 3-dimensional pillar with a square base is discretized into thin sub-volumes along the z axis. We investigate one of such sub-volumes with size ${{L}_{x}}={{L}_{y}}=2000b$ and height $\Delta z=75b$, where b  =  0.256 nm is the norm of the Burgers vector. Initial values for DDD simulations for Systems 1–3 are obtained by randomly distributing ${{N}_{\text{d}}}=50$ circular, discrete dislocation loops with radius R  =  75b inside the volume V (note, that we only chose circular loops for ease of implementation; loops could be arbitrarily shaped—as long as they are closed). Each loop is positioned on a different glide plane and thus loops do not intersect by construction; for the first three systems their centers are distributed on a quadratic domain $0.5~{{L}_{x}}\times 0.5~{{L}_{y}}$ such that all loops are completely contained inside the computational domain. The average density in this volume is ${{\rho}_{0}}=4.8\times {{10}^{15}}$ m⁻². $z$ is the coarse graining height such that the system can be numerically represented as a two-dimensional object—a thin lamella in the x–y plane (figure 1(a)). Continuous field data is obtained by 'smearing-out' the coarse grained fields through the regularization outlined in section 4.2. Initial values for system 4 are obtained by distributing the loop centers on the whole quadratic domain ${{L}_{x}}\times {{L}_{y}}$ while simultaneously taking care of the periodicity. For determining the standard deviation for the Gaussian distribution we estimate the mean dislocation spacing from the density ${{\rho}_{0}}$, which results in a mean dislocation spacing of $s={{L}_{x}}/\left(2\sqrt{2\pi {{N}_{\text{d}}}}\right)\approx {{L}_{x}}/35$. For our examples this estimate of the standard deviation works well for all investigated systems and at all time steps (further details on the choise of s can be found in the appendix A). Only when the mean dislocation spacing changes drastically—in section 5.4 due to dislocation pile ups at grain boundaries—we adjust this value locally during the post-processing step for obtaining reference DDD field data to the actual value of the dislocation spacing. An example for DDD and resulting CDD initial structure that is used for all numerical examples is shown in figure 2: the total dislocation density has a maximum in the center region where most of the loops overlap. At the same time nearly equal numbers of loop segments of all possible orientations are present in the center region which makes the dislocations 'statistically stored'. Hence, the GND density $\mid \boldsymbol{\varrho }\mid$ tends to zero. The curvature density $q$ has the same shape as the total density because initially all loops have the same radius and thus $q=\rho /R$. Although for each of the first 3 studies identical initial values are used, we note that the chosen CDD fields are not special or more suitable than any other set of statistically equivalent initial values.

Zoom In Zoom Out Reset image size

For the numerical solution of the continuum evolution equations a Galerkin finite element scheme together with an implicit time integration scheme is used. To accurately represent derivatives of up to second order, Lagrangian shape functions with cubic polynomials were chosen. The finite element mesh consists of $30\times 30$ elements and thus has an element size which is in the same range as $\Delta z$. When high density gradients at boundaries occur (system 2 and 3) we additionally refine elements in those regions. To increase numerical stability in particular in regions with high fluxes and where the density approaches zero a small viscous term was added to each evolution equation. E.g. in (22) we replace the divergence term $\nabla \cdot \left(v\boldsymbol{\varrho }\right)$ with $\nabla \cdot \left(v\boldsymbol{\varrho }+\epsilon \nabla \boldsymbol{\varrho }\right)$, where ε is chosen sufficiently small such that the physical behavior of the equations is not affected while numerical oscillations are suppressed. The FEM simulations presented below all take less than two minutes computational time on a off-the-shelve computer with a 2.6 GHz dual core processor.

Computation of the evolution of plastic strain can be done for CDD based on the Orowan equation (${{\partial}_{t}}{{\gamma}_{{}}}=\rho bv$) by time integration alongside with the CDD evolution equations (22)–(24). To obtain the initial values for γ from DDD one can compute γ from the GND vector and the relation $\boldsymbol{\varrho }=-\frac{1}{b}\nabla \gamma$: geometrical construction of the swept area for each unit of GND density and subsequent superposition of these areas—followed by division by $\Delta z$—gives the plastic strain. Since the plastic strain is a derived quantity it suffices to show that main CDD field quantities match those from DDD and we do not explicitly show γ.

The DDD model used the representation by chordal splines ($\alpha =1$); for further numerical details on the used Parametric Dislocation Dynamics method please refer to [14, 46].

For this system a constant velocity of v₀  =  0.01 nm/μs is prescribed, and time is expressed as multiple of the duration that density needs to completely traverse the domain, ${{L}_{x}}/{{v}_{0}}$. Boundaries are 'open', i.e. dislocations may leave the volume as if the domain were just a sub-domain of a much larger domain. This case of dislocation loops which have the same curvature and expand with constant velocity is a case for which the CDD theory is supposed to give mathematically exact results, and the system therefore can be used as a reference system for testing the quality of the numerical scheme and if the chosen standard deviation s is appropriate. Initial values are shown in figure 2. Two different simulations are run: discrete microstructure is evolved with the DDD model and—separately—the continuous microstructure is evolved with the CDD model. To compare the two simulations at time t the discrete microstructure was converted into CDD field quantities using the same strategy as for obtaining initial values from discrete data (in the sketch figure 1(f) which is then compared with (e)). Figure 3 shows a snapshot in time, where in the contour plots the dashed lines indicate CDD values and the full lines are the contours of the converted DDD values. We observe that as dislocation loops expand they flow away from the center region where the total density $\rho$ tends to zero (compare $\rho$ in figure 3). The difference in the centers of the loops becomes more insignificant as their radii grow and the resulting density distribution takes a loop-like shape. In each point dislocations now have all approximately the same orientation, they become GNDs which shows in $\rho \approx \mid \boldsymbol{\varrho }\mid$ as can be observed in figure 3. The comparison of DDD and CDD fields shows excellent agreement. Only minor features of the DDD fields are slightly smoothed out which is an artifact of the numerical solution. The volume integrals of $\rho$ and $q/2\pi$ give the total line length and the total number of dislocation loops and are shown in figure 4. For reference we also show the line length for an infinite size system (approximated by a significantly larger domain in x and y direction) for which the total line length can be obtained analytically as $L(t)={{N}_{\text{d}}}2\pi (R+vt)$. Initially, both systems exhibit a linear line length increase (loops expand with constant velocity). As soon as density leaves the volume through the boundaries the line length for the finite size system decreases; the line length in the infinite system, however, continues to increase linearly. CDD data matches the DDD data perfectly and also coincides for the large system with the analytical solution. The right plot of figure 4 shows the total number of dislocation loops. This value stays constant for the infinite system, which is correct since there are neither sources nor annihilation. For the finite system the surfaces act as sinks and the number of loops is reduced. Note that this interplay between density and curvature density will also become relevant when sources and annihilation are to be modeled (which we do not consider within the present work).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Now the outflow of dislocations will be constrained by imposing impenetrable boundary conditions, i.e. the dislocation flux $\rho v$ is enforced to be zero at the boundaries. Numerically, this is done by prescribing a velocity field which decays smoothly to zero directly at the boundary (in physical terms this boundary layer could be imagined as the result of FIB damage in the surface-near regions). To mathematically model the boundary layer we use the following sigmoidal function,

$$\begin{eqnarray}&&f(\xi )=\frac{2}{1+\exp [-\xi \beta ]}-1,\end{eqnarray} \tag{ 27 }$$

where the parameter β controls the shape of the function and is chosen such that the boundary layer is clearly visible and allows to study the dislocation flux in detail ($\beta =20$). ξ is the normalized distance from the boundary, e.g. with $\xi =0\ldots 1$. The velocity field can now be composed from the four contributions

$$\begin{eqnarray}&&v(x,y)={{v}_{0}}\cdot f\left(x/{{L}_{x}}\right)\cdot f\left(1-x/{{L}_{x}}\right)\cdot f\left(\,y/{{L}_{y}}\right)\cdot f\left(1-y/{{L}_{y}}\right).\end{eqnarray} \tag{ 28 }$$

It is depicted in figure 5 and will be used for CDD as well as for DDD simulations. Both simulations are conducted with the initial dislocation structure that was used as well for system 1 (figure 2). Again, the DDD simulation starts from the discrete microstructure, the CDD simulation from the continuous microstructure as obtained by D2C conversion. Both simulations evolve independently. Only in a post-processing step we again convert the resulting DDD microstructure for reference purposes into continuous field quantities as outlined before. A snapshot in time of the respective evolved dislocation field quantities is shown figure 6 with CDD values shown as dashed contour lines and converted DDD field values shown as solid lines (also see the figure in the appendix A and the supplementary movie 1 (stacks.iop.org/MSMSE/23/085003/mmedia) which additionally shows the evolution of average curvature $k=q/\rho$). During the simulation dislocation loops initially expand freely until parts of them reach the boundary layer within which the velocity decays to zero. The non-zero velocity gradient rotates line segments such that they eventually are parallel to the boundary, i.e. their line tangents are aligned with the contours of the velocity field. This causes the initial SSD density to increasingly become geometrically necessary when dislocations pile up against the boundaries (compare $\rho$ and $\mid \boldsymbol{\varrho }\mid$ in figure 6). Additionally, dislocations need to bend with a high curvature radius near the corners of the domain. This shows in the high value of the curvature density in these regions. For a more detailed view figure 7 shows horizontal and diagonal line plots through the domain for two different time steps. We observe again an excellent agreement between the DDD and CDD field values for all times, and only small deviations can be seen close to the boundaries and in particular for the curvature density.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The total line length inside the system serves as a good plausibility test: initially we again would expect free loop expansion until finally all loops are aligned with the boundary and are (nearly) rectangular in shape. Figure 8 shows the total line length versus time, and it is obvious that at all times DDD and CDD agree very well: both approach the horizontal tangent when all lines are aligned with the boundaries. The fact that DDD values are slightly below the horizontal tangent is correct because loops are not exactly rectangular but still have small rounded corners. For the DDD curve we directly used the line length from the discrete data while for the CDD data we first had to integrate $\rho$. The small difference in the final value of the line length arises from discretization errors due to very steep gradients at the boundaries. In this respect, more realistic systems that also consider elastic dislocation interactions will numerically behave even better because pile ups always have a finite size.

Zoom In Zoom Out Reset image size

The third system studies the idealization of a plastically deforming matrix into which an elastic, circular inclusion is embedded. Dislocations can move freely in the matrix while the inclusion acts as an impenetrable obstacle for dislocation motion. Annihilation of lines is not considered, hence, the number of dislocation loops stays constant; for real systems this is a strong restriction, but again, the philosophy is to decouple interactions and reactions from the purely kinematic effects. The computational model is—as in system 1—a quadratic domain (${{L}_{x}}={{L}_{y}}$) with open boundaries and contains the inclusion in the center. The inclusion itself is represented by an internal boundary which we model as a circular sub-domain (radius R  =  L_x/10) with zero velocity in the center region and a smooth transition from v  =  0 to a constant velocity v₀ away from the inclusion (see figure 9(a)). As a result of the velocity field dislocations are first slowed down on their way towards the inclusion and eventually freeze in their motion so that they cannot enter the inclusion. Both simulation methods use the same initial values as before. This does place some dislocations inside the inclusion which is unphysical. It turned out, however, that the influence on the resulting microstructure is only small. Therefore, we decided to consistently use the same initial values as for all other systems as well. Evolving the DDD and CDD system gives typical dislocation patterns as shown in figure 10, additional plots including error plots can be found in the appendix A.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

An auspicious feature is the formation of pile-ups of bent lines around the inclusion: when expanding dislocation loops touch the inclusion they get a 'dent' (see the colored line in figure 10(a)) and adjust to the shape of the inclusion. At the same time, line segments further away start 'spiraling' around the inclusion and deposit an increasing number of segments around the inclusion (also see the figure in the appendix A and the supplementary movie) (stacks.iop.org/MSMSE/23/085003/mmedia). This behavior has a number of consequences and can be observed in the CDD representation as well:

• Circular pile-ups of dislocation density can be seen in the upper middle and right plots of figure 10 : $\rho$ has a high value in the vicinity of the inclusion because density from piled-up dislocations together with pre-existing line segments in the inside of the precipitate add up. In conjunction with the GND density $\mid \boldsymbol{\varrho }\mid$ we infer that dislocations in the inside are mainly statistically stored while pile-ups around the inclusion are—literally—geometrically necessary.
• Circular GNDs around the inclusion form an inverted loop, see figure 10 bottom right and middle plot. Taking a look at the sign of e.g. the edge GNDs ${{\varrho}_{\text{e}}}$ around the inclusion it is found that their orientation was rotated by 180° (also compare the initial values in figure 2). For instance the upper half is now negative and therefore would move downwards but is hindered by the inclusion.
• Bent dislocations around the inclusion have negative curvature—the above statement about an inverted loop already suggests that the curvature of dislocations around the inclusion should be negative which in fact can be observed in the two plots on the right of figure 11 where the average curvature $k=q/\rho$ was computed and plotted along two lines passing through the inclusion. The inclusion has a radius of L_x/10 which is equivalent to a curvature of ${{k}_{\text{inc}}}=0.020$ nm⁻¹. Due to the smooth interface the velocity is zero only at about L_x/16 which is equivalent to a curvature of ${{k}_{\text{min}}}=0.031$ nm⁻¹. In the bottom right plot of figure 11 we find the negative curvature peak values at $\mid {{k}_{\text{CDD}}}\mid \approx 0.040$ nm⁻¹ in good agreement with both the geometry value ${{k}_{\text{min}}}$ and the DDD curvature value of $\mid {{k}_{\text{DDD}}}\mid \approx 0.028$ nm⁻¹.

Zoom In Zoom Out Reset image size

Deviations between DDD and CDD are larger than for the previous two benchmark models. They are particularly large in those regions which exhibit large gradients in the velocity and in the CDD field values, both of which cause increasing discretization errors. What is the average response of this system? Comparing the temporal evolution of the DDD line length in figure 10 to system 1 shows that the inclusion creates additional density which is clearly visible from approximately t  =  0.15 ${{L}_{x}}/{{v}_{0}}$ on and is caused by lines bending around the inclusion. CDD is lagging behind and is initially even slightly below the line length of system 1 (due to the lower average velocity as compared to system 1). After $t\geqslant 0.25$ ${{L}_{x}}/{{v}_{0}}$ the CDD model is then also showing an additional increase in line length.

The difference between CDD and DDD is the result of larger errors close to the inclusion where CDD almost always underestimates the true density values. Increasing the resolution of the spatial discretization of the finite element scheme does not yield any appreciable alleviation. The reason is the very complex deformation state of dislocations that brings the CDD theory to its limits because near the inclusion the line curvature can be different for different line orientations in the same averaging volume—a detail that cannot be represented with this simplified variant of CDD (this CDD theory only can represent one average curvature value for each point/averaging volume). Although a number of details of this system have been predicted properly—e.g. the line curvature around the inclusion which is important when it comes to hardening effects in terms of line tension—this particular CDD formulation only roughly predicts the line length production close to the inclusion. More elaborate evolution equations which contain more information about the orientation distribution and curvature of dislocations are required as e.g. those derived in [45, 50]. The D2C strategy for more detailed CDD formulations, however, will still remain the same.

The previous benchmark systems have in common that they start with random initial values which during time evolution become strongly polarized through the geometrical constraints imposed by boundary conditions (e.g. dislocations adjust their line orientation and curvature to the shape of external or internal boundaries and become geometrically necessary). System 4 will now investigate how good CDD performs in a bulk-like situation where the total density consists to a large extend of statistically stored dislocation density superimposed with (smaller) GND density fluctuations. Periodic boundary conditions are used to eliminate geometrical constraints. Initial values consist of a statistically homogeneous random distribution of 200 dislocation loops of the same radius, and their centers are distributed across the whole volume ${{L}_{x}}\times {{L}_{y}}\times \Delta z$. This results in the same average density ${{\rho}_{0}}$ as in the center region of Systems 1–3; all other parameters stay the same.

The velocity and loops' curvature radius are constant everywhere. Since this study analyzes the quality of the representation of flow and evolution of density fluctuations, it is useful to decompose the total density into fluctuating and volume-averaged contributions, $\rho (x,y)={{\rho}_{\text{fluct}}}(x,y)+\langle \rho \rangle$. For plotting the field data we additionally normalize the fluctuations with the variance of the respective field, e.g. in figure 12 we plot for the total density $(\rho -\langle \rho \rangle )/\text{var}(\rho )$. Figure 13 shows line plots of initial and evolved DDD and CDD data, additional error plots can be found in the appendix A. The supplementary movie S3 shows the time evolution of density fields. Clearly, all details of the discrete microstructure are accurately reproduced. Fluctuations evolve in a non-trivial manner and we attribute negligible deviations to numerical errors. The average behavior matches exactly the DDD values. The overall very good agreement is not surprising after seeing that already system 1, which initially had a similar random structure in the center of the domain, gave very good agreement between DDD and CDD. It is, however, important to show that this good agreement is sustained during time evolution and that even the increase in the range of the fluctuations are properly represented. An application based on this type of system has been studied in the context of dislocation patterning [41] where in particular the fluctuations together with the correct line length increase and a Taylor-type flow stress were identified as key ingredients for pattern formation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size