Simulations of grain growth in tungsten armor materials under ARC plasma edge operation conditions using an integrated plasma-edge/materials model

TL;DR

This paper presents a new integrated model to simulate grain growth in deuterium-exposed tungsten under plasma conditions, revealing rapid grain growth near plasma and stability inside the material, with implications for tungsten's use in fusion reactors.

Contribution

The study introduces a two-dimensional vertex dynamics model fitted to atomistic data for simulating grain growth in tungsten under plasma exposure, a novel approach for this application.

Findings

Fast grain growth occurs near plasma at 1400 K within 100 seconds.

Microstructure remains stable deep inside the material at 1000 K over several days.

Conventional tungsten is highly susceptible to grain growth above 1000 K.

Abstract

An integrated model of grain growth deuterium-exposed tungsten polycrystals, consisting of a two-dimensional vertex dynamics model fitted to atomistic data, has been developed to assess the grain growth kinetics of deuterium-exposed polycrystalline tungsten (W). The model tracks the motion of grain boundaries under the effect of driving forces stemming from grain boundary curvature and differential deuterium concentration accumulation. We apply the model to experimentally-synthesized W polycrystals under deuterium saturated conditions consistent with those of the ARC concept design, and find fast grain growth kinetics in the material region adjacent to the plasma (at 1400 K, <100 seconds for full transformation), while the microstructure is stable deep inside the material (several days to complete at a temperature of 1000 K). Our simulations suggest that monolithic W fabricated using conventional techniques will be highly susceptible to grain growth in the presence of any driving force at temperatures above 1000 K.

Figures

Tables

Table 1. Table 1: Parameters used to define the mobility tensor in W for different misorientations along the [110] tilt axis. The superscript ‘P’ denotes pristine W, while ‘H’ refers to hydrogen-loaded GB [ 39 ] .

$\theta$ [deg]	$f_0^{\text{P}}$ [s-1]	$Q_a^{\text{P}}$ [eV nm-1]	$f_0^{\text{H}}$ [s-1]	$Q_a^{\text{H}}$ [eV nm-1]	$h$ [Å]	$b$ [Å]	$A$ [nm2]
26.52	$4.0\times {10}^{17}$	1.99	$5.6\times {10}^{14}$	1.68	1.58	0.8	29.20
58.99	$4.9\times {10}^{17}$	11.09	$8.4\times {10}^{15}$	10.26	2.60	2.3	27.17
65.47	$1.0\times {10}^{16}$	6.87	$3.3\times {10}^{12}$	6.31	2.54	2.0	30.92
70.53	$1.7\times {10}^{14}$	1.95	$1.3\times {10}^{13}$	1.56	1.30	0.9	28.94
77.88	$2.0\times {10}^{15}$	2.58	$1.0\times {10}^{14}$	2.55	2.66	1.5	42.62
82.95	$1.8\times {10}^{16}$	1.00	$3.0\times {10}^{15}$	1.00	2.66	1.3	25.25
86.62	$2.3\times {10}^{18}$	1.11	$3.1\times {10}^{16}$	1.18	2.66	1.0	29.32

Equations

$$\vec{v}=\mathbf{M}\vec{f}=\left(\begin{array}{ccc}{M}_{11}& M_{12}& M_{13}\\ M_{12}& M_{22}& M_{23}\\ M_{13}& M_{23}& M_{33}\end{array}\right)\left(\begin{array}{c}\psi \\ {\sigma }_{12}\\ {\sigma }_{13}\end{array}\right)$$

$$v_1=M_{11}\psi +M_{12}{\sigma }_{12}$$

$$v_2=M_{21}\psi +M_{22}{\sigma }_{12}$$

$$M_{11}=M(\theta ,T)n_1$$

$$M_{21}=M(\theta ,T)n_2$$

$$v_1=M_{11}\psi$$

$${\vec{f}}_{\gamma }=-\frac{\gamma (\theta )}{R}\vec{n}$$

$$E_i^{\mathrm{tot}}=AL\left(c_1+c_2\right)\Delta h_{\text{D}}$$

$$E_f^{\mathrm{tot}}=A\left(c_1(L+\Delta x)+c_2(L-\Delta x)\right)\Delta h_{\text{D}}$$

$$f_{\text{D}}^{\text{tot}}=-\frac{\Delta E^{\text{tot}}}{\Delta x}=-\frac{E_{\text{f}}^{\text{tot}}-E_{\text{i}}^{\text{tot}}}{\Delta x}=A\left(c_1-c_2\right)\Delta h_{\text{D}}$$

$${\vec{f}}_{\text{D}}=\Delta h_{\text{D}}(c_1-c_2)\vec{n}=\Delta h_{\text{D}}\Delta c\vec{n}$$

$${\vec{f}}_{\psi }={\vec{f}}_{\gamma }+{\vec{f}}_{\text{D}}=\left[-\frac{\gamma (\theta )}{R}+\Delta h_{\text{D}}\Delta c\right]\vec{n}.$$

$${\vec{v}}_{\mathrm{TJ}}={\mathcal{M}}_{\mathrm{TJ}}\sum _{k=1}^3{\gamma }_k{\vec{t}}_k$$

$$\delta V=\frac{V_{\text{GB}}-N{\Omega }_a}{2A}$$

$$c_i=\frac{n_{\mathrm{D}}}{V_i}$$

$$V_f=V_i\pm A\delta V$$

$$c_f=\frac{n_{\mathrm{D}}}{V_f}=\frac{n_{\mathrm{D}}}{V_i+A\delta V}=\frac{n_{\mathrm{D}}}{A\left(w\pm \delta V\right)}$$

$$\Delta c=c_f-c_i$$

$$=c_i\frac{\delta V}{w\pm \delta V}\approx \pm c_i\frac{\delta V}{w}$$

$$f_{\text{RSW}}(x,a)=\sin \left(\frac{\pi }{2x}\right)\left[1-a\log \left\{\sin \left(\frac{\pi }{2x}\right)\right\}\right]$$

$$\mathbf{M}=\frac{f_{0}A}{kT}\sum _m\left({\vec{H}}^m\otimes {\vec{H}}^m\right)\exp \left(-\frac{Q_a^m}{kT}\right)$$

$$f_A=1-\exp \left(-a(t-t_0{)}^n\right)$$

Full text

Keywords: Tungsten; Fusion materials; Plasma-facing materials; Grain growth; Modeling and simulation; Multiscale modeling

1 Introduction

Structural materials for fusion energy systems must be designed to face extreme conditions during reactor operation [1, 2, 3, 4]. In particular, divertors in magnetic fusion devices will be subjected to nominal heat fluxes exceeding 10 MW$\cdot$m*-2* [5, 6, 7, 8] (and even higher during sporadic discharge events known as edge localized modes, which can reach heat loads of up to 100 MW$\cdot$m*-2* [9, 10, 11]). Tungsten (W) has emerged as the main candidate plasma-facing component (PFC) material due to a set of favorable properties, including a high melting point (3695 K), high thermal conductivity (173 W$\cdot$m*-1*$\cdot$K*-1* at room temperature), low tritium storage, low coefficient of thermal expansion ($4.3\times 10^{-6}$ K*-1*), low neutron activation, low erosion rate and high-temperature yield strength. However, the use of W is limited by a narrow operational temperature window, with its upper limit set by recrystallization, and its lower one by a generally unacceptably-high brittle-to-ductile (DBTT) transition temperature [12, 13] (which can range between room temperature for high-purity, cold-rolled specimens [14, 15, 16] to over 500∘C for conventionally forged specimens [17, 18, 19]). In addition, due to the intermittent –but cyclic– nature of plasma kinetic behavior during fusion reactor operation, it is also important to consider elements of thermal fatigue in the mechanical performance of W components [20, 21, 22, 23, 24, 25, 26, 27].

The DBTT of conventional polycrystalline tungsten depends on the material microstructure, and can be lowered by plastic deformation to values even below room temperature [28, 14, 29]. Pre-straining by, e.g., cold working or rolling is the basis for the processing of the so-called ‘fusion-grade’ tungsten [30, 31]. However, the deformation structures resulting from low-temperature deformation are thermally unstable, and predispose the material to a set of restoration processes during operation at high temperatures: recovery, recrystallization, and grain coarsening. The unavoidable occurrence of one or more of these processes during fusion device operation ultimately brings W back to its inherently brittle undeformed state, effectively setting practical limits on the material’s lifetime. Recrystallization and grain growth are enabled by grain boundary (GB) motion, driven by differential driving forces resulting from chemical potential jumps across the boundary [32], or from the anisotropy of the GB energy itself [33]. A natural strategy to delay grain coarsening processes is then to enable Zener pinning by selective alloying and GB segregation [34, 35] (notwithstanding other detrimental processes that might result from alloying). For example, the recrystallization temperature limit in W-3Re alloys doped with K was seen to increase beyond 850∘C relative to that of pure W [35].

Under operation, the divertor will be inevitably exposed to hydrogen isotopes (H/D/T) from the plasma, resulting in the penetration of H atoms inside the material and interactions with its internal microstructure. In particular, grain boundaries are known to trap H ions with energies ranging between 0.8 and 1.2 eV [36, 37, 38], a process which has been seen to decrease their intrinsic mobility [39, 40] and alter the equilibrium concentration profiles that would be expected from Fickian kinetics. However, the collective effect of trapped and free H isotopes on microstructural evolution during grain growth is not well understood. Past studies [41, 42, 43, 44] have shown no clear difference between pristine and H-loaded W specimens at temperatures near the recrystallization limit. Still, characterization of internal H concentration profiles during H-plasma expoure is exceedingly challenging, and the quantitative correlation between H-isotope concentration and grain growth remains to be established. It is thus crucial to gain a detailed understanding of the processes leading to microstructural changes of W-based components exposed to hydrogen isotopes.

In this paper, we carry out vertex dynamics simulations of grain growth in plasma-exposed polycrystalline W by adapting an existing model for recrystallization [32] to the particular conditions expected in Commonwealth Fusion Systems’ next-of-a-kind fusion concept ARC [45]. The ARC is in an early design stage, so the choice of the surface and substrate material, as well as the details of the geometric design and allowable plasma parameters are still under development. These simulations are thus intended as a guide to setting reasonable limits and expectations on materials performance under ARC’s planned operational conditions. As boundary conditions, we use depth-resolved temperature distributions and D-T concentrations recently calculated by an integrated plasma-materials interaction model. The model captures boundary plasma conditions as obtained by the UEDGE code[46], ion-surface interactions by F-TRIDYN, and sub-surface gas dynamics by Xolotl [47]. This integrated approach was then used to simulate the exposure in ten locations along the outer divertor target in detached plasma conditions [48].

The paper is organized as follows. First, a review of the theory behind grain boundary evolution models and the implementation adopted in this work are provided in Sec. 2. The results of the paper are presented in Sec. 3, including the selection of the microstructure from characterization maps and grain simulations under several different scenarios. We provide a discussion of our findings in Sec. 4, and finish with the main conclusions of the paper in Sec. 5.

2 Theory and methods

2.1 Vertex model for grain boundary kinetics

Our model of grain boundary evolution is based on a two-dimensional (2D) vertex model of the polycrystal microstructure [32]. The model consists of physical nodes, representing triple junctions (TJ), and virtual nodes connected as a piecewise discretization of a grain boundary. A summary of the features of the model is provided in the next section.

2.1.1 Virtual node evolution equations:

The motion of a grain boundary (GB) is defined by a generalized displacement $\mathbf{u}=(u_{1},u_{2},u_{3})$ with component $u_{1}$ representing motion in the direction of the GB normal, and $u_{2}$ and $u_{3}$ representing in-plane motion. The GB velocity, $\vec{v}=\dot{\vec{u}}$, responds to a set of driving forces defined by the vector $\vec{f}$:

[TABLE]

where $\bm{M}$ is the mobility matrix. The first component of $\vec{f}$, $f_{1}=\psi$, includes all ‘chemical’ sources of motion. The second and third components of $\vec{f}$ are $f_{2}=\sigma_{12}$ and $f_{2}=\sigma_{13}$, which are shear stresses along the $x_{2}$ and $x_{3}$ directions in the GB plane. These stresses lead to so-called coupled shear boundary motion through the corresponding mobility components $M_{12}$, $M_{23}$, and $M_{13}$.

In 2D, eq. (1) is reduced to:

[TABLE]

For zero stress conditions (see Secs. 2.1.3 and 2.1.4) and assuming a decomposition of the GB normal $\vec{n}$ as $n_{1}=\cos\phi$ and $n_{2}=\sin\phi$ (where $\phi$ is the polar angle formed by the GB normal relative to a 2D cartesian frame of reference), the corresponding mobilities can be expressed as:

[TABLE]

where $M_{11}$ captures motion of a GB vertex along the GB normal, and $M_{21}$ represents its motion in the GB plane. $M(\theta,T)$ is defined in Sec. 2.3.3. Motion governed by $M_{11}$ responds to a mechanical driving force, while motion defined by $M_{21}$ reacts to ‘configurational’ driving forces. In this work, we disregard the latter, which leads to the following simplified equation of motion:

[TABLE]

2.1.2 Chemical driving force for GB motion:

The components of $f_{1}=\psi$ in eq. (1) that contribute to GB motion include:

• (i)

Curvature-driven motion:

[TABLE]

where $R$ is the radius of curvature of the approximate circle representing the GB at a given instant in 2D, and $\gamma(\theta)$ is the misorientation-dependent GB energy ($\theta$ is the misorientation angle), $\vec{n}$ denotes the unit normal vector of the GB surface, oriented outward in the direction of the local curvature.

• (ii)

Differential deuterium concentration: If we consider the ideal bi-crystal shown in Fig. 1, composed of two identical grains of length $L$, a GB area $A$, and each containing a deuterium concentration $c_{1}$ and $c_{2}$, a force to move the GB towards the grain with the higher deuterium concentration will exist. This force can be derived by calculating the total energy of the deuterium in solution in the bi-crystal:

[TABLE]

where $\Delta h_{\text{D}}$ is the enthalpy of solution (formation enthalpy) of deuterium in tungsten per atom, taken here as constant. After the GB moves by a distance $\Delta x$ into grain 1, the current total energy of bi-crystal is

[TABLE]

For infinitesimal displacements of the GB, $c_{1}$ and $c_{2}$ can be assumed to remain unchanged after the GB moves, the driving force for the GB motion resulting from deuterium concentration jump is equal to:

[TABLE]

The above expression is for a true force acting on a point. Dividing this force by the GB area $A$ gives the capillary pressure due to the concentration differential across the interface:

[TABLE]

where $\Delta c$ is the deuterium concentration jump across the GB and $\vec{n}$ is the unit normal vector of the GB. A positive force is taken for $c_{1}>c_{2}$, i.e., with the GB moving in the direction of grain with the higher deuterium concentration.

Combining both forces above, the total chemical driving force $\vec{f}_{\psi}$ experienced by a GB is:

[TABLE]

which is used in eq. (2) in its scalar form $\psi=\|\vec{f}_{\psi}\|$. The function $\gamma(\theta)$ is defined in Sec. 2.3.1. For $\Delta h_{\text{D}}$, we take a value of 1.4 eV consistent with DFT calculations [49, 50].

2.1.3 Shear-coupled motion:

The shear-coupling factor (the ratio of the shear velocity to the GB migration velocity) has two components $\beta_{k}=v_{k}/v_{1}$ corresponding to orthogonal shear directions contained in the GB facet. This factor can be obtained by measuring $v_{1}$ and $v_{k}$ under an applied shear stress $\sigma_{1k}$, i.e., $\beta_{k}=M_{kk}/M_{1k}$, which can be done directly using MD simulations [51, 52, 53, 54].

In the present study no stress is applied and thus there is no direct shear coupled GB motion from external stresses. Rather, it is intrinsic (i.e., ‘regular’) GB motion that induces internal shear stresses. As discussed in Sec. 2.1.4, these stresses transfer to triple junctions, creating the possibility for twin or dislocation nucleation into the grain [55]. However, at present we assume that triple junctions have sufficiently high mobilities to relieve these stresses as they appear.

2.1.4 Triple junction motion:

The role of triple junctions (TJ) in controlling grain growth in annealed polycrystalline materials was recognized in early theoretical and simulation studies [56, 57, 58]. Specifically, when grain growth in a polycrystal is governed by the mobility of triple junctions, the kinetics deviates from the classical Johnson-Mehl-Avrami [59] –or its variant, the Mullins-Neumman [58]– evolution law.

Triple junctions move according to the following viscous law:

[TABLE]

where the index $k$ refers to the three GB that converge on the TJ and ${\cal M}_{\rm TJ}$ is the triple junction mobility. $\gamma_{k}$ is the energy per unit area of GB $k$. The three forces from the GBs converging into the TJ are oriented along the GB tangents $\vec{t}_{k}$. In equilibrium, $\vec{v}_{\rm TJ}$ is zero and if the GB energy is independent of misorientation, one gets the dihedral angle condition (Herring relation). Indeed, eq. (9) implies that grains with dihedral angles of less than $90^{\circ}$ are always unstable and will eventually disappear.

When the ratio $\Lambda=L{\cal M}_{\rm TJ}/M_{\rm GB}$ is large ($L$ is the grain size: note that the units of the TJ and GB mobilities differ by one length dimension), the dihedral angle condition is maintained because the TJ can adjust instantaneously to the motion of the GBs. When $\Lambda$ is small, steady state dihedral angles that deviate from the static solution may be seen, but those are physical and may ultimately control microstructural evolution. In other words, if triple junctions can move sufficiently rapidly to remain in equilibrium with respect to the GBs, then the microstructure evolves by pure chemical forces with fixed angle junctions. This type of motion may be viewed as a two-step process whereby (i) GBs migrate in accordance with eq. (1), thereby pulling the TJs out of equilibrium, and (ii) the TJs migrate to restore the equilibrium angles following eq. (9). Experimental observations and simulation show that TJs do indeed have finite mobility, as indicated by steady-state dihedral angles of migrating TJs that differ from their equilibrium values [60, 57, 56, 55].

Quantitative assumptions to define the value of ${\cal M}_{\rm TJ}$ and $\Lambda$ are provided in Sec. 2.3.4.

2.2 Deuterium concentration jump across a grain boundary

The presence of grain boundaries is known to alter impurity concentration profiles relative to those in a homogeneous (e.g., single crystal) environment. Studies have revealed discontinuous segregation patterns across twin boundaries [61], thin films [62], and grain boundary complexions [63]. When the solute concentration is uniform, segregation enrichment or depletion are common in polycrystals [64, 65]. However, when the concentration field is biased (e.g., subjected to a gradient), ‘jumps’ in solute concentration can be seen across boundaries without changing the concentration gradient. These jumps lead to differential solute concentrations that result in a driving force for GB motion. The process is schematically shown in Figure 2, where a solute concentration profile is superimposed on a generic polycrystalline system, leading to jumps in concentration across each GB. In this section we derive the relationship between the solute concentration differential, $\Delta c$, and the GB misorientation, $\theta$.

Our approach is based on the consideration of an excess ‘free’ volume, $\delta V$, of a GB relative to the bulk, an idea going back to solute segregation models by Aaron and Bolling in 1972 [66]. In the absence of more accurate ways to obtain it, Wolf proposed a direct correlation between $\delta V$ and the GB misorientation following the Shockley-Read expression for GB energies [67]. However, the excess volume can now be readily obtained from atomistic simulations as:

[TABLE]

where $V_{\text{GB}}$ is the volume of a simulation box containing a periodic bi-crystal (i.e., two equivalent grain boundaries with equal surface area $A$), $N$ is the total number of atoms in the simulation, and $\Omega_{a}$ is the atomic volume (equal to $a_{0}^{3}/2$ in BCC metals, where $a_{0}$ is the lattice constant). Equation (10) compares the actual volume of a system containing a fixed number of atoms and two grain boundaries with the volume that those atoms would occupy in a bulk system. To remove the scale effect associated with the arbitrary value of $N$, $\delta V$ is normalized by the GB area (and thus expressed in units of length). Note that the excess volume can be positive or negative, depending on where the tilt angle is measured from.

The concentration jump across a GB with an associated excess volume of $\delta V$ can be obtained by assuming that the impurity concentration in the bulk experiences a In the case of deuterium atoms, the same number of particles, $n_{\rm D}$, is now distributed over a different volume than in the bulk. The concentration prior to encountering a GB, $c_{i}$, can be written as:

[TABLE]

where $V_{i}=wA$ is the volume of a region adjacent to the GB with a width equal to $w$. At the GB, the total available volume changes to:

[TABLE]

where $A=L^{2}$ is the GB surface area. The deuterium concentration after the GB is now:

[TABLE]

The concentration jump, $\Delta c=c_{f}-c_{i}$, is then:

[TABLE]

where we have assumed that $\delta V\ll w$. Thus, the excess volume leads to a correction in the concentration profile that scales with the concentration measured at the GB on the side closest to the exposed wall of the material.

2.3 Physical parameters

2.3.1 Grain boundary energies:

In this paper, we take GB energies with a non-monotonic dependence on misorientation obtained using atomistic calculations. Specifically, we consider energies for $[110]$ symmetric tilt boundaries calculated using semi-empirical potentials [68, 69, 70] validated with electronic structure calculations [71, 72]. The $\gamma_{[110]}$-$\theta$ function for these family of GB is presented in Figure 3a.

The data displayed in the figure are then used to fit a set of Read-Shockley-Wolf functions (RSW) functions [74]:

[TABLE]

where $x=\frac{\theta-\theta_{\text{min}}}{\theta_{\text{max}}-\theta_{\text{min}}}$, $\theta$ and $a$ are the misorientation angle and shape parameter [75] respectively, while $\theta_{\text{min}}$ and $\theta_{\text{max}}$ represent the periodicity range of the $\gamma_{[110]}(\theta)$ function (0 to $90^{\circ}$ in this case). These two parameters are obtained by fitting eq. (12) to atomistic calculations.

In our two-dimensional continuum grain growth model, the misorientation is calculated as the angular offset between the projected crystal axes of the grains on both sides of the GB. The inclination is then taken as the angel between the GB plane normal, $\vec{n}$, and these axes.

2.3.2 Grain boundary excess volume:

The excess volume for symmetric tilt boundaries in W has been calculated using atomistic methods by Zheng et al. [73]. The values for $[110]$-tilt GB are given in Figure 3b normalized by the value of $a_{0}$ in W (3.15 Å). In agreement with the Wolf model [67], the dependence of $\delta V$ on $\theta$ closely mimics that of $\gamma_{\rm[110]}$ with $\theta$ (Fig. 3a). As Fig. 3b shows, this dependence is non-monotonic and here we simply perform a linear interpolation between data points for a given value of $\theta$. The data also show that $\delta V$ is only a fraction of $a_{0}$, ranging between zero and 0.14.

Several studies have estimated the width of grain boundaries in metals to be in the 0.5-to-1.0-nm range [76, 77]. For that reason, the parameter $w$ in eq. (11) is set to $2a_{0}$ here. This choice of $w$ justifies the approximation $\delta V\ll w$ adopted in Sec. 2.2.

2.3.3 Grain boundary mobilities:

Grain boundary mobility is one of the most challenging aspects of polycrystal evolution. While it is reasonable to assume that bi-crystallography is the primary determinant of grain boundary mobility (through its five-dimensional relationship [60]), recent mobility calculations obtained by matching polycrystal evolution models to dynamic grain boundary evolution observations have yielded no definitive correlation between misorientation, inclination, and GB motion [78]. This is consistent with MD studies in Ni bicrystals [79], while in others, such as Al [80, 81], mobilities display a marked misorientation dependence. In the present study, the problem is compounded by the presence of deuterium, which may hinder grain boundary mobility relative to the pristine material condition.

To address these issues, here we use recent atomistic results of GB mobility in W [39]. In particular, we use a definition of $\bm{M}$ based on the decomposition of a general grain boundary into a set of $m$ disconnections described by a structure vector $\vec{H}$ [53]:

[TABLE]

where $f_{0}$ is a prefactor, $k$ is Boltzmann’s constant, $\vec{H}=[b,h,0]$ ($b$ and $h$ are the magnitude of the disconnection’s Burgers vector and step height respectively), $Q_{a}$ is the activation energy, and $A$ is again the GB area. In the current model, we take only the leading term of the sum in eq. (13), i.e., $m=1$. The parameters obtained by molecular dynamics simulations are listed in Table. 1, where results for pure (pristine) W are given as well as for hydrogen-loaded GBs with H concentration to be 0.01 Å*-2* (labeled with superindices ‘P’ and ‘H’, respectively).

Surface plots of $M_{11}$ for the pristine and hydrogen-loaded cases as a function of misorientation and temperature are given in Fig. 4. As the figures show, GB mobility is virtually zero at temperatures below 900 K and for misorientations below 75∘. In the range where the mobilities are non-negligible, deuterium-loaded GB display values that are 10 times lower than those calculated for pristine W.

2.3.4 Triple junction mobility:

Here we consider two regimes for TJ motion:

The first regime is characterized by a ratio $\Lambda<1$, generically representing low TJ mobility. Accordingly, the kinetics is expected to be controlled by TJ –not GB– motion, with stresses arising from shear-coupled motion building up at unrelaxed TJs. As indicated above, beyond a certain stress threshold, such buildup may result in nucleation of twins or dislocations from the TJ into the grain. However, at present, this effect is not captured in the simulations. 2. 2.

The second scenario involves $\Lambda>1$, representing conditions when GBs are expected to control the evolution of the microstructure. In this case, shear-coupled stresses are immediately relaxed by TJ motion, resulting in no emission of twins or dislocations.

In the results that will be presented below, we parametrically vary $\Lambda$ between 0.1 and 100, and study its effect on the resulting microstructures. For clarity, we use a value of $L$ equal to one in units such that $\Lambda$ reflects the non-dimensional ratio of mobilities in a direct way.

3 Results

3.1 Deuterium concentration profile and temperature distribution.

Figure 5 shows depth profiles for the deuterium concentration and the temperature of the outer divertor target after exposure to 10 cycles of 900 s plasma with a 60-s break, as calculated by Xolotl for the location of peak plasma and heat fluxes at several locations near the strike point. Figure 5a reveals a flat steady-state deuterium concentration profile beyond the first 0.1 microns of material, at a level of $\approx 1.9\times 10^{24}$ m*-3*. Thus, in our simulations we simply consider a uniform D distribution throughout the entire depth of the material (subject to the concentration jumps across GB discussed in Sec. 2.2). For its part, the temperature follows a decreasing linear profile from which a constant gradient of $3.2\times 10^{5}$ K$\cdot$m*-1* can be extracted. In our simulations, we fix the temperature at the edge of the computational domain to 1370 K, and apply the constant gradient to set the temperature across it.

3.2 Material microstructure and selection of the computational domain

Next, we create synthetic two-dimensional microstructures from EBSD (electron backscattering diffraction) characterizations of pure W specimens obtained by the following sequence: (1) cold isostatic pressing at room temperature; (2) sintering at 1800$\sim$2200∘C; (3) hot rolling at 1400$\sim$1600∘C to an area reduction of approximately 80%; (4) stress relief heat treatment at 900∘C for 20 min. Figure 6a shows the ‘as-fabricated’ microstructure colored according to the grain orientation as indicated in the attached standard triangle color map. The figure shows the surface exposed to D from the plasma (l.h.s. of the image) and a scale marker for reference. Figure 6b shows the color map in the stereographic triangle used to represent the crystal orientation in each grain. Based on the color map, the material appears to be textured with a slight predominance of crystal orientations near the [101] crystallographic axis (colored in red in the image).

The entire microstructure shown in Fig. 6a contains 2994 grains and over 18,000 grain boundaries, which is presently beyond the computational capabilities of our vertex dynamics model. Instead, we use a subset of the configuration shown in the figure for the actual simulations. The actual microstructure employed in the simulations is highlighted by a white dashed box at the bottom left corner of Fig. 6a (referred to as ‘inset 1’). This substructure contains approximately 115 grains and 584 grain boundaries, which makes it tractable computationally. However, to ensure that this subset is representative of the full set we conduct a statistical verification exercise comprising the grain size distribution and the misorientation distribution function (MDF), i.e., the normalized misorientation angle distribution. Figure 6c shows a histogram with the initial grain size distribution of the microstructure shown in Fig. 6a. The figure also shows a log-normal fit to the histogram data, from which an average grain size of 5.0$\sim$7.0 microns can be extracted. Importantly, the log-normal fit to the grain distribution extracted from the inset is also shown in the figure and found to be in reasonable agreement with that corresponding to the entire microstucture.

The MDF is shown in Figure 6d. The MDF displays a bimodal shape, with a peak at misorientations of $\theta\approx 10^{\circ}$ and another at $\approx 45^{\circ}$. As in the case of the grain size distribution, the MDF is fitted to a dual log-normal composite function, corresponding to the entire and subset structures shown in Fig. 6a, are in good agreement with one another. Combined, the results in Figs. 6c and 6d provide evidence of the statistical equivalence between the full microstructure and the subsets chosen for computational efficiency.

A similar analysis performed for the substructure highlighted by the white dashed box in the bottom right portion of the EBSD image (to be used for simulations described below in Sec. 3.3.3) led to the same conclusions.

Figure 7a shows the grain substructure of subset ‘1’ (also colored according to the map provided in Fig. 6b), while Fig. 7b shows the corresponding grain boundary misorientation map. As the subfigure shows, the microstructure in subset ‘1’ is dominated by low-angle grain boundaries ($<20^{\circ}$). Next, we carry out simulations of microstructure evolution in this subset with the fully-parameterized model.

3.3 Grain growth simulations

3.3.1 Microstructure evolution:

Here, we track the temporal evolution of the system and compute the fraction of transformed surface area by discretizing the entire grain domain into a $100\times 100$ uniform grid. A given location within the grid is said to be transformed if a GB sweeps over it at any point in time. Each grid point is counted only once as transformed, regardless of whether it has been swept more than once. The fraction of the swept area, $f_{A}(t)$, is then calculated as the ratio of transformed grid points relative to the total number.

Figure 8 shows a sequence of snapshots of the microstructure in subset ‘1’ from $t=0$ until the maximum transformation is achieved. Visually, the smaller grains in the lower third sector of the image are seen to undergo rapid growth, relative to the rest of the substructure. As well, grains with crystal orientations close to the $[101]$ zone axis (colored red in the figure) appear to suffer a slower transformation. The evolution of $f_{A}(t)$ with time is shown in Figure 8f, showing saturation after 100 s of simulated time. The simulations predict that grain growth will begin instantaneously upon exposure of W to deuterium at the temperatures of the first wall (Fig. 8a). Rapid growth persists during the following 10 s of evolution (Figs. 8b-c), primarily by consumption of small grains by larger grains everywhere in the specimen. Subsequently, the system experiences slower growth (Figs. 8d-e), likely due to a ‘locked’ microstructure limited by the slow growth of only one or very few specific grains. This manifests itself as a ‘bend’ in the growth curve between 40 and 70 s in Fig. 8f. The system is seen to saturate at approximately 100 s of evolution, with 80% of the total area swept.

3.3.2 Quantitative analysis of the grain substructure

The general observations from the system evolution presented in the prior section are that the microstructure coarsens by absorption of smaller grains by larger ones. Here we provide a more quantitative analysis of the main microstructural descriptors and their evolution during the grain growth process.

Figure 9 contains information about the evolution of the grain size with time. The grain size distribution at three stages during the growth process, 0, 60, and 90% completed, is shown in Figure 9a. Clearly, the distributions become more uniform, deviating from the original log-normal. This is accompanied by coarsening, as shown in Fig. 9b, where a monotonic reduction in the total number of grains from the original 115 to merely 41 of them. Likewise, the average grain size increases from 10.5 to $\approx 18$ microns.

The change in the misorientation distribution function between $f_{A}=0$ and $f_{A}=0.9$ is given in Figure 10. As the figure shows, while the changes are not drastic, a depletion in the relative proportion of LAGB in favor of higher angle ones can be clearly observed. A priori, this may suggest that highly mobile GB are activated towards eliminating grains with a high proportion of LAGB (which are less mobile according to Fig. 4). To elucidate this point, however, a more appropriate analysis is to directly quantify the relative participation of each GB in the entire grain growth process. To that effect, in Figure 11a we show the distribution of misorientations at the beginning of the simulations ($F_{2}$, equivalent to Fig. 6c) next to their contribution to the swept area fraction, $F_{1}$.

The relative weight of each misorientation bin in the total grain growth process is then established by obtaining the ratio between $F_{1}$ and $F_{2}$, which is shown in Figure 11b. It is clear that HAGB, $>50^{\circ}$, contribute the most proportionally to grain growth, almost a factor of two more than any other GB misorientation. This points to GB mobility as the controlling factor.

3.3.3 Low temperature grain growth simulations:

Last, we carry out simulations of grain growth at 1000 K for the microstructure contained in inset ‘2’ shown in Fig. 6a. Such microstructure satisfies the statistical equivalence conditions of grain size and grain boundary misorientation distributions as described in Sec. 3.2, but we use it instead of inset ‘1’ for statistical variation. The results are shown in Figure 12, where the fraction of swept area is shown as a function of time. The microstructure begins to evolve after 3 s of clock time, but its transformation is not completed until after $10^{6}$ s (approximately 11.5 days). The curve for the microstructure in inset ‘1’ (adjacent to the plasma) is also shown for comparison, together with the case for pristine W where only capillary forces (GB curvature) are at play. Clearly, temperature is the most significant factor for microstructural evolution and grain boundary activation.

4 Discussion

4.1 Model assumptions and simplifications

• •

Triple junction mobility: Here, we have assumed that TJ mobility is much higher than GB mobility such that triple junctions immediately adopt their equilibrium configuration in response to GB motion. TJ mobility is a challenging parameter to obtain experimentally, as triple-junctions are of course 2D representations of 3D structures [56]. There are very few modeling works [82, 83]. Recent works have established the effect of TJ mobility on grain growth (described by a growth scaling law of \sim$$t^{n}) when considered in conjunction with other parameters [84, 85]. In particular, these are expected values of $n$ depnding on what controls grain growth kinetics:

• (i)

Controlled only by GB mobility: a time exponent of $n=1/2$ is to be expected.

• (ii)

Controlled by the energy of the grain boundaries together with the mobility of the triple junctions: $n=1$.

• (iii)

Controlled by the energy of the triple junctions together with the mobility of the grain boundaries: $n=1/3$.

• (iv)

Finally, when the kinetics is solely controlled by the triple junction mobility, $n=1/2$ again.

Of course, other factors such as defects, dislocations, impurities, etc, lead to deviations from these expected values. In our study, our model runs under conditions (i) and (ii) above, which is consistent with the Avrami exponents that will be discussed in Sec. 4.2.

• •

Two-dimensional model: The present model considers polycrystals as 2D systems representing an axis zone projection of full 3D microstructures. The model is suitable for reconstructions of polycrystals from EBSD maps provided that the Euler angles for all grains and the zone axis are known (both of which are the case here). The implications of representing 3D polycrystals in 2D have been discussed in our previous work [32]. Suffice it to say that Potts Monte Carlo simulations in 2D and 3D have shown that $n=1/2$ (representative of 2D kinetics) can be considered to be the asymptotic long-time kinetic growth exponent even for three-dimensional grain growth [86].

• •

Selection of $[110]$ tilt axis: The reason to consider GB properties for $[110]$ tilt boundaries in this study is simply that that orientation is the only one for which GB energies and mobilities are jointly available. It is not clear at this point what effect, if any, this has on grain growth and microstructural evolution. Given that GB energies for $[110]$ tilt boundaries are quantitatively similar to other crystal axes (see, e.g., ref. [68]), and that the dependence of GB mobilities on temperature and misorientation are not expected to differ substantially from those used here (see ref. [39]), we speculate that the choice made here is reaosnably representative of all GB tilt orientations. It has been shown, however, that the structure of any GB with the misorientation axes $[110]$ (and $[111]$) may be represented by a limited number of basic structural units [87], such that $[110]$ tilt boundaries are an ideal choice for general grain growth simulations.

• •

Influence of plasticity and residual stresses: As indicated in Sec. 1, fusion-grade tungsten typically involves warm or hot rolling to bring the DBTT to engineering acceptable levels. This will inevitably lead to the buildup of dislocation density and the onset of residual stresses. In our previous work [32], we have shown how to capture the effect of plasticity in driving grain growth processes, and, per eq. (1), the model is intrinsically capable of incorporating stresses. However, the microstructure maps featured in Fig. 6a do not include any information beyond grain orientations and sizes. Thus, the present work does not intend to be a numerical analysis of the response of realistic tungsten to fusion conditions, but a study focused on the effects of deuterium exposure to pure W polycrystals.

Despite the simplification and assumptions above, the model has value in assessing the susceptibility of monolithic W to grain growth under ARC operational conditions, as well as estimating the time scales over which it takes place, as will be shown below.

4.2 Avrami exponents and dimensionality of the kinetics

To assess the kinetics of grain growth in an analytical fashion, next we fit the curves in Fig. 10 to a modified Johnson-Mehl-Avrami-Kolmogorov (JMAK) model [88, 89], i.e., where the fraction of transformed area can be represented by the following expression:

[TABLE]

where $a$ is a fitting constant, $n$ is the the Avrami exponent (which describes the dimensionality of the governing mechanisms behind grain growth), and $t_{0}$ represents the incubation time before the start of the transformation. For a grain growth process purely defined by GB motion, a value of $n\approx 2$ is expected. As indicated in Fig. 10, the Avrami exponents extracted from our simulations are 0.73 for the cases containing deuterium and 1.13 for the case without. These values are in excellent agreement with experimentally-measured values. Alfonso et al. [90] found a value of $n=1.1$ in their study of recrystallization kinetics of warm-rolled tungsten, while Wang et al. [91] found values between 0.78 and 1.0 in heavily pre-deformed W at intermediate temperatures. Finally, Shah et al. [41] found an average value of 1.4 in their study of H-plasma exposed W specimens. Again, this is remarkable agreement with our simulation results given the differences and assumptions pertaining to the model.

Furthermore, tests with varying temperatures allow for the determination of the activation energy for grain growth, which, for H-loaded samples, was estimated to be approximately 4.5 eV [41]. Referring to this value to the activation energies GB listed in Table 1 suggests that the evolution of the system is controlled by GB with misorientation angles of $\approx 78^{\circ}$. This is consistent with the results presented in Fig. 11b, where HAGB clearly contribute the most to the entire grain growth process.

4.3 Effect of deuterium on grain growth

As implemented in the model, D has two main effects on microstructural evolution. The first one is preferential segregation towards grain boundaries, which results in concentration jumps that scale with misorientation according to eq. (11) and Fig. 3b. This leads to differential D accumulation across GB, which, in turn, results in driving forces for GB evolution. The second is a reduction in GB mobility due to an increased activation energy as indicated in Table 1.

From the simulations, it is clear that D has a delaying effect on GB kinetics, as manifested by the lower Avrami exponent of 0.72 compared to a value ¿1.0 for pristine W. Interestingly, there is no effect on the incubation time (i.e., the time at which the transformation commences), which is apparently controlled only by temperature.

4.4 Implications for plasma-facing armor

Ultimately, our simulations suggest that monolithic tungsten fabricated using conventional techniques will be highly susceptible to grain growth in the presence of any driving force at temperatures above 1000∘C. Here, we have considered GB curvature (probably established during grain deformation processes due to hot rolling) and differential deuterium accumulation at GB (assumed to occur due to excess atomic volume in the vicinity of boundaries) as the main driving forces in our study. We have shown that those alone are sufficient to drive microstructural evolution at 1400 K to ¿80% transformed fraction within less than 100 s. Decreasing the temperature to 1000 K, delays the onset of grain growth by four orders of magnitude in time, clearly highlighting the critical role of temperature in driving internal changes.

However, one must also consider that, in general, pre-strained tungsten will contain a relatively high density of residual dislocations. Since plastic deformation is highly-dependent on crystal orientation, it is expected that large differential dislocation density accumulation around GB will result in large driving forces for microstructural evolution, leading to even faster transformations. A possible technique to delay grain growth and stabilize pre-deformed W microstrutures is to take advantage of Zener pinning by alloying with selected solutes or seeding with second-phase particles. In fusion-grade W, promising systems include W-Re [92, 93, 94, 95], W-K [96], W-Re-K [35, 97], or W infiltrated with lanthanide oxide particles [98, 99]. In most of these cases, it is indeed found that W alloys display a higher recrystallization temperature relative to pure W.

Finally, under irradiation, two competing effects can be envisaged to appear. First, an extra driving force is introduced in the form of irradiation defects accumulating in different proportions on both sides of the GB. This could presumably accelerate microstructural evolution in the same manner as GB curvature, plastic deformation, and deuterium exposure do. At the same time, however, irradiation defects are known to strongly interact with grain boundaries, leading to complex structures [100] which may pin GB and hinder the onset of grain growth or alter the GB mobility[101]. From the limited experimental evidence currently available [102, 35, 103], it appears that irradiation damage also leads to delayed microstructural evolution processes in comparison with unirradiated W. However, we leave the effects of irradiation on grain growth for future studies.

5 Conclusions

In summary, we have developed a vertex dynamics model in 2D to assess the grain growth kinetics of deuterium-exposed polycrystalline tungsten. The model tracks the motion of grain boundaries under the effect of driving forces stemming from grain boundary curvature and deuterium concentration jumps. All the model parameters have been obtained from atomistic simulations, including grain boundary energies, mobilities, excess volumes, and deuterium solution energies. The simulations assume saturation conditions for D-exposed tungsten as obtained from integrated plasma (UEDGE code) -material (F-TRIDYN and Xolotl) modeling of the outer divertor target exposed to 10 cycles of 900 s of plasma exposure with 60-s breaks, at locations of peak plasma and heat fluxes near the strike point.

Our main findings are listed below:

• •

The initial material microstructure has been analyzed statistically from two-dimensional EBSD maps. We find that the grain size distribution follows a log-normal distribution with an average grain size of approximately 11 microns. The grain boundary misorientation function is bimodal, with one peak at 10∘ and another at 45∘.

• •

In the plasma-facing region (first 100 microns, 1400 K), grain growth commences immediately upon deuterium saturation of the microstructure, taking less than 100 s to complete. The average grain size is seen to grow from 11 to 18 microns.

• •

We find that high-angle grain boundaries ($\approx 60^{\circ}$) contribute proportionally the most to grain growth.

• •

Avrami exponents of microstructural evolution laws range between 0.72 for the case of D-exposed tungsten and 1.2 for unexposed W. These are in very good agreement with experimental measurements.

• •

Temperature has a drastic impact on grain growth kinetics, delaying the incubation time from $1.0$ to $10^{4}$ s when the temperature is reduced from 1400 to 1000 K.

• •

Deuterium has a pinning effect on grain boundaries, reducing their mobility and delaying the onset and progression of grain growth.

Acknowledgments

This work was funded under the INFUSE program, a DOE-FES public-private partnership, under contract no. DE-SC0024662 between the University of California Los Angeles, the University of Tennessee-Knoxville, and Commonwealth Fusion Systems. The authors also acknowledge support from DOE-FES under contracts DE-SC0018410 (J.Y. and J.M.) and DE-SC0023180 (S.B., A.L., and B.D.W.). N.M. gratefully acknowledge support from the U.S. DOE, Office of Science, Office of Fusion Energy Sciences, and Office of Advanced Scientific Computing Research through the Scientific Discovery through Advanced Computing (SciDAC) project on Plasma-Surface Interactions (Award No. DESC0008875). We thank Shuhei Nogami from ALMT for providing the tungsten material and Blake Emad at University of North Texas for performing the EBSD analysis on the tungsten specimens. This manuscript has been assigned LA-UR-25-27324.