A multi-scale model for stresses, strains and swelling of reactor components under irradiation

One of the challenges associated with the design, construction and operation of a fusion or an advanced fission power plant is the need to predict stresses, strains, and swelling resulting from the exposure of structural components of a reactor to irradiation during operation [1–5]. These stresses and strains have a microscopic origin and stem from the fact that radiation defects have substantial elastic relaxation volumes [6, 7], which give rise to strong local deformation of the lattice. For example, the elastic relaxation volume of a self-interstitial atom defect in tungsten, predicted by density functional theory, is $\Omega_{\rm rel}=1.67\Omega_0$, where $\Omega_0$ is the volume of an atom, whereas the relaxation volume of a vacancy is $\Omega_{\rm rel}=-0.37\Omega_0$, see e.g. [7, 8]. The large positive mismatch between elastic relaxation volumes of self-interstitial and vacancy type defects in metals, which accumulate in the microstructure as a result of irradiation, gives rise to the local volumetric expansion, and produces stresses in materials exposed to irradiation, resulting in macroscopic swelling and heterogeneous deformation of reactor components.

Models for swelling developed since the 1970s focus on one particular aspect of the phenomenon, the diffusion-mediated preferential absorption of self-interstitial atom defects by dislocations [6, 9–11]. There is extensive literature on the dynamics of accumulation of defects in microstructure and particularly on the evaluation of dislocation bias factors [6, 12–17], which are then used as input parameters for mean-field rate theory models describing the dynamics of growth of voids in materials exposed to irradiation [9–11, 18]. These mean-field models do not require, and do not evaluate, the macroscopic stresses and strains resulting from the accumulation of defects. Neither do they take into account the fact that elastic relaxation of complex defect configurations involves the accumulation of very strong local deformation of the lattice near the defects, with a consequence that the relaxation volume of a complex defect cluster is not equal to the sum of relaxation volumes of constituting defects. The defect configurations formed in collision cascades initiated by high-energy neutrons [19–24] have a fairly complex structure, and their relaxation volumes, and local strains and stresses are sensitive to the relative proximity of individual defects.

A particularly significant practical aspect of the problem of accumulation of radiation defects in structural materials exposed to high-energy neutron irradiation during power plant operation, which has so far eluded a satisfactory solution, is related to the fact that in a treatment of deformation of a component exposed to irradiation, not only the presence of defects in the bulk of the material, but also boundary conditions at surfaces, define the final configuration of stresses and strains. For example, it is known that the magnitude of the total elastic relaxation volume of a point defect depends on boundary conditions at surfaces of the material where it resides and, surprisingly, the role of boundary conditions does not diminish in the macroscopic limit [25, 26].

Below we develop a real-space multi-scale model for strains and stresses produced in reactor components by spatially distributed complex configurations of radiation defects produced by irradiation. We show that if the properties of materials are reasonably isotropic, it is the defect relaxation volume density, i.e. the real-space distribution of relaxation volumes of self-interstitial and vacancy defects and clusters of such defects in a reactor component, that fully defines the resulting macroscopic strains and stresses. Individual defects and clusters of defects act as sources of strains and stresses, which are determined self-consistently, taking into account the appropriate boundary conditions. The fundamental notions relating macroscopic strains and stresses to microstructure are the dipole tensor and relaxation volume of a defect object, which can be an individual point defect or a large cluster of such defects formed in a collision cascade. We show that the dipole tensor of a defect, which has so far only been applied to nano-scale point defects [7, 25, 27], can in fact be used for characterizing defect objects of arbitrary size. This is a mere consequence of the fact that there is no characteristic spatial scale in elasticity equations. At a large distance from an arbitrarily complex configuration of defects, for example, a cluster of defects containing thousands of individual point defects, the elastic strain generated by the cluster equals

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon_{ij}({\bf r})=-P_{kl}{\partial ^2\over \partial x_j\partial x_l}G_{ik}({\bf r}-{\bf R}), \label{elastic_strain} \nonumber \end{align} \tag{ 1 }$$

where ${\bf R}$ is the position of the defect cluster, $G_{ik}({\bf r}-{\bf R})$ is Green's function of elasticity equations [28], and P_kl is the dipole tensor of the defect configuration. Equation (1) fully defines the strain field generated by a defect cluster at distances that are much greater than its characteristic spatial extent.

It is well established how to compute elastic dipole tensors of point defects [7, 27, 29]. Below, we show how elastic dipole tensors can be defined and evaluated for defect clusters produced by the collapse of entire high-energy collision cascades. The condition that makes it possible to apply the dipole tensor concept to the treatment of strains and stresses in a reactor component is that the spatial scale of a defect object in the radiation-induced microstructure is always small in comparison with the spatial scale of a component.

Historically, elastic dipole tensors were introduced in order to treat long-range elastic fields of nano-scale point defects, and elastic interactions between such defects, as well as between defects and dislocations [25, 30, 31]. Still, the fact that the universal scale-free nature of elasticity equations enables extending the notion of elastic dipole tensor of a defect to mesoscopic and macroscopic scales, does not appear to have been recognized.

In the treatment of a spatially-localized defect configuration, which may include defects produced by the collapse of an entire cascade or even a group of cascades, it is only necessary to take into account the discreteness of the atomic lattice in the most strongly distorted core regions of defect structures. At large distances from a defect configuration, lattice discreteness is no longer significant and the rate of variation of atomic displacements as a function of spatial coordinates is small $|\partial u_i/\partial x_j|\ll 1$. In this limit, the strain and stress fields generated by a defect configuration are fully characterized by a symmetric $3\times 3$ tensor P_kl. The volume dipole tensor density, and the density of relaxation volumes of defects, vary from one point in the material to another, depending on the local degree of exposure of a material to irradiation. The dipole tensor can be defined independently of the nature of a defect object, for example, it can be defined for a point defect, a dislocation loop [7, 27, 31], or a fairly complex configuration of atomic displacements produced by a collision cascade as a whole. In the latter case, although the volume of the spatial region involved in the calculation of P_kl contains a large number of atoms, it is still small in comparison with the size of a reactor component, justifying the fundamental approximation on which equation (1) is based.

In the next section we outline the formalism involved in the evaluation of stress and strain fields associated with spatially distributed defect structures. We derive the fundamental equations defining strains, stresses and displacements in an arbitrary volume of a material containing spatially distributed defects. If the material is elastically isotropic and defects have no preferred orientation, the strains and stresses are fully determined by the distribution of spatially-varying density of elastic relaxation volumes, a dimensionless function that acts as a source term in equations for the irradiation-induced stresses and strains. Then we show how the relaxation volumes of individual point defects can be computed using ab initio density functional theory methods, and provide accurate values of relaxation volumes computed for several metals with cubic crystal structure. Subsequently, we compare accuracy of ab initio and atomic relaxation methods using semi-empirical interatomic potentials, and show how to obtain accurate values of dipole tensors for entire configurations of defects produced by the collapse of collision cascades initiated by high-energy neutrons. For both cases, we discuss a practical way of defining the density of relaxation volumes of defects and relate it to a measure of exposure of a material to neutron irradiation. We then give a comprehensive derivation of the macroscopic elasticity equations including the density of body forces associated with the accumulation of defects, and provide a full exact solution of these equations for the case of a spherical shell with a source of irradiation at its centre. The solution makes it possible to determine swelling, strains and stresses everywhere in the shell, and it also illustrates the pivotal part played by the boundary conditions for the elasticity equations. Finally, we discuss a finite element implementation of the method and examine a possible application of our approach to the in silico evaluation of operational performance of a fusion power plant.

At a large distance from a strongly deformed atomic configuration, for example a defect or a cluster of defects, the field of displacements of atoms from their equilibrium positions in the lattice is given by the equation [25]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u_i({\bf r})=-P_{kl}{\partial\over \partial x_l}G_{ik}({\bf r}-{\bf R}).\label{displacements} \nonumber \end{align} \tag{ 2 }$$

In this equation, P_kl is the elastic dipole tensor of the defect configuration that, according to (2), fully defines the elastic field that the configuration generates in the material. $G_{ik}({\bf r}-{\bf R})$ is the elastic Green's function [28], which in the isotropic elasticity approximation has the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle G_{ik}({\bf r})={1\over 16\pi \mu (1-\nu)r} \left[(3-4\nu)\delta _{ik}+ {x_ix_k\over r^2}\right].\label{elasticGF} \nonumber \end{align} \tag{ 3 }$$

Here, μ is the shear modulus of the material and ν is the Poisson ratio, see [32]. From (2) and (3) we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u_i({\bf r})&={1\over 16\pi \mu (1-\nu)|{\bf r}-{\bf R}|^2} \nonumber \\ &\times\left[2(1-2\nu)P_{il}\eta_l-\eta_iP_{kk}+3\eta_i(\eta_k P_{kl} \eta_l)\right],\label{displacements1} \nonumber \end{align} \tag{ 4 }$$

where summation over repeated indexes is assumed, and $\newcommand{\e}{{\rm e}} \eta_i=({\bf r}-{\bf R})_i/|{\bf r}-{\bf R}|$.

Equations similar to (4) have been extensively explored in connection with the treatment of elastic fields of point defects on the nano-scale [7, 27, 33, 34]. However, as the form of equation (4) suggests, there is no specific spatial scale at which it should be applied, as elastic fields have no intrinsic spatial scale associated with them. For example, using formula (4) one can approximate the elastic field of a fairly large defect structure just as easily as one can do this for a point defect.

The magnitude and the angular anisotropy of elastic field is described by the matrix elements of tensor P_kl, which can be evaluated using atomistic simulations with periodic boundary conditions, using the equation [7, 25, 29]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{kl}=-\int_V \sigma_{kl}({\bf r}){\rm d}^3r=-V{\overline \sigma}_{kl},\label{evaluationPij} \nonumber \end{align} \tag{ 5 }$$

where $V$ is the volume of the simulation cell containing the defect structure, and ${\overline \sigma}_{kl}$ is the volume average stress tensor. In the literature [25, 29], equation (5) has so far been only applied to point defects, and hence the size of the cell used in simulations was relatively small, in most cases not exceeding 10³ atoms. However, since the notion of the dipole tensor defined above is entirely general, equation (5) can be applied to a defect structure of arbitrary size, potentially involving millions of atoms and many individual point defects and clusters of defects.

It is convenient to represent P_kl in the following form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{kl}=C_{klmn}\Omega_{mn}, \label{dipoleTensorFromStrain} \nonumber \end{align} \tag{ 6 }$$

where $\Omega _{mn}$ is an auxiliary $3\times 3$ tensor, also uniquely characterizing the defect structure, and C_ijkl is the tensor of elastic constants. The advantage that this representation offers is that it provides a simple and robust way of computing the relaxation volume of an arbitrary complex defect configuration as a whole. Indeed, the relaxation volume of a defect configuration in a body with surfaces free of tractions equals [25]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{relaxationVolume} \Omega_{\rm rel}=S_{mmkl}P_{kl}= \Omega_{mm} ={\rm Tr}\,\Omega, \nonumber \end{align} \tag{ 7 }$$

where $S_{ijkl}=(C_{ijkl}){}^{-1}$ is the elastic compliance tensor [35]. Tensor $\Omega_{mn}$ can be expressed in terms of its eigenvectors and eigenvalues as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{partialRelaxationVolume} \Omega _{mn}=\sum_{s=1}^3 \Omega^{(s)}e^{(s)}_me^{(s)}_n, \nonumber \end{align} \tag{ 8 }$$

where eigenvectors ${\bf e}^{(s)}$ form a mutually orthogonal set, and eigenvalues $\Omega ^{(s)}$ have the meaning of partial relaxation volumes. They satisfy the condition

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Omega_{\rm rel}=\sum_{s=1}^3\Omega^{(s)}, \nonumber \end{align} \tag{ 9 }$$

which, if the defect configuration has been formed as a result of a collision cascade event initiated by an energetic neutron, gives the total relaxation volume of cascade debris. Swelling, often occurring under neutron irradiation, results from the accumulation of uncompensated relaxation volumes of defects and defect clusters produced in collision cascades initiated by neutrons.

In engineering applications, where the scale of effective volumes containing defects is large, and defects and defect clusters adopt random spatial orientations, the eigenvectors of tensor $\Omega_{mn}$ may point in any spatial direction with equal probability. For example, consider the relaxation volume tensor of a crowdion defect [7] shown schematically in figure 1

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Omega_{mn}=\Omega^{(1)}l_m l_n+{\Omega^{(2)}\over 3}\delta_{mn}. \nonumber \end{align} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

The direction vector ${\bf l}$ of the axis of the defect can adopt any of the four $\langle 1\, 1\, 1\rangle$ type directions. As a result, the average over defect orientations relaxation volume tensor is diagonal with respect to m and n.

Similarly, we can define an average over orientations for the relaxation volume tensor of a complex defect

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle \Omega _{mn}\rangle =\sum_{s=1}^3 \Omega^{(s)}\langle e^{(s)}_me^{(s)}_n\rangle.\label{average_orientation} \nonumber \end{align} \tag{ 11 }$$

Here $\langle ...\rangle$ denotes averaging over all the directions in the solid angle of $4\pi$, namely

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \langle\,f(\theta, \phi) \rangle ={1\over 4\pi}\int f(\theta, \phi) \sin \theta {\rm d}\phi {\rm d}\theta. \nonumber \end{align*}$$

For a unit vector in 3D space ${\bf e}=(\sin \theta \cos \phi, \sin \theta \sin \phi, \cos \theta),$ we find that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle e_i\rangle &=0, \nonumber \\ \langle e_i e_j\rangle &= {1\over 3}\delta _{ij}, \nonumber \end{align} \tag{ 12 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle \Omega_{mn}\rangle ={1\over 3} \Omega_{\rm rel}\delta_{mn}. \nonumber \end{align} \tag{ 13 }$$

Substituting this into equation (6) we arrive at

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle P_{kl}\rangle ={\Omega_{\rm rel}\over 3}C_{klmm}, \nonumber \end{align} \tag{ 14 }$$

which for crystals of cubic symmetry becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle P_{kl}\rangle ={\Omega_{\rm rel}\over 3}(C_{11}+2C_{12})\delta_{kl}. \nonumber \end{align} \tag{ 15 }$$

In the isotropic elasticity limit the above equation for the dipole tensor reduces to [7, 31]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle P_{kl}\rangle ={2\mu \over 3}\Omega_{\rm rel}{1+\nu \over 1-2\nu}\delta_{kl}=B\Omega_{\rm rel}\delta_{kl},\label{isotropic} \nonumber \end{align} \tag{ 16 }$$

where B is the bulk modulus of the material.

In applications, where macroscopic elastic fields are formed as a result of superposition of microscopic fields created by millions of individual defects and defect clusters, which all adopt arbitrary orientations in the lattice, it is often sufficient to use the above isotropic approximation for the dipole tensor. The advantage offered by this approximation is that a defect object is characterized by only one parameter, its relaxation volume $\Omega_{\rm rel}$. In some cases, for example in crystals with non-cubic symmetry or in the presence of significant external anisotropic elastic strain field, it may prove necessary to take into account the anisotropy of defect structures.

The field of elastic displacements generated by a defect object, the dipole tensor of which has the form (16), equals

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u_i({\bf r})={\Omega_{\rm rel}\over 12\pi}\left({1+\nu \over 1-\nu}\right) {({\bf r}-{\bf R})_i\over |{\bf r}-{\bf R}|^3}. \nonumber \end{align} \tag{ 17 }$$

The strain field, corresponding to this field of displacements, is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon_{ij}({\bf r})={\Omega_{\rm rel}\over 4\pi |{\bf r}-{\bf R}|^3}{1+\nu \over 1-\nu}\left({1\over 3}\delta_{ij}-\eta_i\eta_j\right).\label{distortion2} \nonumber \end{align} \tag{ 18 }$$

To compute the elastic stress associated with the strain field given by (18), we need to multiply the strain tensor by the four-index tensor of elastic constants, which in the isotropic elasticity approximation equals [26]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C_{ijkl} =\mu {2\nu \over 1-2\nu} \delta_{ij} \delta_{kl} + \mu \left ( \delta_{ik} \delta_{jl}+\delta_{il}\delta_{jk} \right).\label{elastic_constants_isotropic} \nonumber \end{align} \tag{ 19 }$$

The resulting expression for the stress field generated by a defect object at ${\bf R}$ is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{ij}({\bf r})={\mu \Omega_{\rm rel}\over 2\pi |{\bf r}-{\bf R}|^3}{1+\nu \over 1-\nu}\left({1\over 3}\delta_{ij}-\eta_i\eta_j\right), \; {\bf r}\ne {\bf R}.\label{stress2} \nonumber \end{align} \tag{ 20 }$$

We now define the central notion of our treatment, the density of relaxation volumes of defects, which is a dimensionless quantity equal to the total relaxation volume of all the defect configurations in a unit volume of material

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \omega _{\rm rel}({\bf r})=\sum _{a}\Omega _{\rm rel}^{(a)}\delta ({\bf r}-{\bf R}_{\alpha}),\label{relaxation_volume_density} \nonumber \end{align} \tag{ 21 }$$

where summation is performed over all the defect clusters a in the material. Because of the properties of the Dirac delta-function, if we integrate $\omega_{\rm rel}({\bf r})$ over a certain volume element in the material, only the defect clusters contained in that volume element are going to contribute to the integral. There is a fundamental relation between the notion of the number of displacements per atom [36], which is a measure of exposure of a material to irradiation, and the density of relaxation volumes of defects that such exposure produces. Below, we show that function $\omega _{\rm rel}({\bf r})$ provides a quantitative measure of the magnitude of sources of stresses and strains produced in materials by radiation defects.

A quantity similar to (21) was introduced by Caturla et al [37, 38] in the context of modelling microstructural evolution and swelling of materials under irradiation. To explain swelling, in [37, 38] the density of relaxation volumes of vacancies was taken as a positive quantity. This does not agree with ab initio calculations, showing that vacancies have negative relaxation volumes whereas self-interstitial atom defects have positive relaxation volumes. Swelling results from the agglomeration of self-interstitial defects into dislocation loops and extended dislocation network; vacancies aggregate into voids, the relaxation volume of which in the macroscopic limit is negative and small [39].

The elastic strain, generated by all the defects produced by irradiation and distributed with certain density in the bulk of a reactor component, can be computed by integrating the point source function (18) with the density of relaxation volumes of defects (21) as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon_{ij}({\bf r})={1\over 4\pi}{1+\nu \over 1-\nu}\int {\omega_{\rm rel}({\bf R}) \over |{\bf r}-{\bf R}|^3}\left({1\over 3}\delta_{ij}-\eta_i\eta_j\right){\rm d}^3R.\label{strain_integrated} \nonumber \end{align} \tag{ 22 }$$

We note that equation (18) follows from (22) if we assume $\omega_{\rm rel}({\bf R})=\Omega_{\rm rel}\delta({\bf R}-{\bf R}_0)$.

Computing the stress generated in a material by a continuous distribution of defects involves an element of mathematical subtlety. First, we re-write equation (22) in the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon_{ij}({\bf r})=-{1\over 12\pi}{1+\nu \over 1-\nu}\int \omega_{\rm rel}({\bf R}){\partial ^2\over \partial x_i \partial x_j} {1\over |{\bf r}-{\bf R}|}{\rm d}^3R.\label{strain_integrated1} \nonumber \end{align} \tag{ 23 }$$

To find the elements of the stress tensor

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{ij}({\bf r})=C_{ijkl}\epsilon_{kl}({\bf r}), \nonumber \end{align} \tag{ 24 }$$

we multiply (23) by the four-index tensor of elastic constants (19). This gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{ij}({\bf r})=&-{\mu\over 6\pi}{1+\nu \over 1-\nu}\int \omega_{\rm rel}({\bf R}) {\partial ^2\over \partial x_i\partial x_j}{1\over |{\bf r}-{\bf R}|}{\rm d}^3R\nonumber \\ &-{\mu \over 12\pi}\left({1+\nu \over 1-\nu}\right)\left({2\nu \over 1-2\nu}\right)\delta_{ij}\nonumber \\ &\times \sum_k \int \omega_{\rm rel}({\bf R}) {\partial ^2\over \partial x^2_k}{1\over |{\bf r}-{\bf R}|}{\rm d}^3R, \label{stress_integrated1} \nonumber \end{align} \tag{ 25 }$$

where the second term, proportional to $\delta_{ij}$, arises from the first term in (19). Noting that

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \sum_k {\partial ^2\over \partial x^2_k}{1\over |{\bf r}-{\bf R}|}=\Delta {1\over |{\bf r}-{\bf R}|}=-4\pi \delta({\bf r}-{\bf R}), \nonumber \end{align*}$$

we find the stress produced by a continuous distribution of defects

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{ij}({\bf r}) =&-{\mu\over 6\pi}{1+\nu \over 1-\nu}\int \omega_{\rm rel}({\bf R}) {\partial ^2\over \partial x_i\partial x_j}{1\over |{\bf r}-{\bf R}|}{\rm d}^3R.\nonumber \\ &+{\mu \over 3}\left({1+\nu \over 1-\nu}\right)\left({2\nu \over 1-2\nu}\right)\delta_{ij}\omega_{\rm rel}({\bf r}). \label{stress_integrated2} \nonumber \end{align} \tag{ 26 }$$

The density of body forces arising from a continuous distribution of defects can now be found by differentiating the stress tensor with respect to x_j, namely [40]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_i({\bf r}) &= -{\partial \sigma _{ij}({\bf r})\over \partial x_j} ={2\mu \over 3}\left({1+\nu \over 1-\nu} \right){\partial \over \partial x_i}\omega_{\rm rel}({\bf r})\nonumber \\ &-{\mu \over 3}\left({1+\nu \over 1-\nu} \right)\left({2\nu \over 1-2\nu} \right){\partial \over \partial x_i}\omega_{\rm rel}({\bf r}). \label{body_force} \nonumber \end{align} \tag{ 27 }$$

Deriving the first term in the above equation involved integration by parts

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &-{\mu\over 6\pi}{1+\nu \over 1-\nu}\int \omega_{\rm rel}({\bf R}) {\partial \over \partial x_i}{\partial^2\over \partial x_j^2}{1\over |{\bf r}-{\bf R}|}{\rm d}^3R\nonumber \\ &={\mu\over 6\pi}{1+\nu \over 1-\nu}\int \omega_{\rm rel}({\bf R}) {\partial \over \partial x_i}[4\pi \delta ({\bf r}-{\bf R})]{\rm d}^3R\nonumber \\ &=-{\mu\over 6\pi}{1+\nu \over 1-\nu}\int \omega_{\rm rel}({\bf R}) {\partial \over \partial X_i}[4\pi \delta ({\bf r}-{\bf R})]{\rm d}^3R\nonumber \\ &={2\mu\over 3}{1+\nu \over 1-\nu}{\partial \over \partial x_i} \omega_{\rm rel}({\bf r}). \label{stress_differential_equation} \nonumber \end{align} \tag{ 28 }$$

The addition of the two terms in (27) gives the density of body forces resulting from the accumulation of defects in the material

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_i({\bf r})=-{2\mu\over 3}\left({1+\nu \over 1-2\nu}\right){\partial \over \partial x_i}\omega_{\rm rel}({\bf r}).\label{body_force1} \nonumber \end{align} \tag{ 29 }$$

This equation is one of the key findings of this work.

For completeness, we also note the expression for the field of elastic displacements produced by the defects distributed in the bulk of a reactor component with relaxation volume density (21), namely

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u_i({\bf r})&=-{1\over 12\pi}\left({1+\nu \over 1-\nu}\right){\partial \over \partial x_i}\int {\omega _{\rm rel}({\bf R}) \over |{\bf r}-{\bf R}|} {\rm d}^3R\nonumber \\ &={1\over 12\pi}{1+\nu \over 1-\nu}\int \omega _{\rm rel}({\bf R}) {({\bf r}-{\bf R})_i\over |{\bf r}-{\bf R}|^3} {\rm d}^3R.\label{displacements_integrated} \nonumber \end{align} \tag{ 30 }$$

Equations (22)–(25) fully define the elastic field generated directly by the sources of stresses and strains associated with relaxation volumes of spatially distributed radiation defects. However, these equations do not take into account the boundary conditions at surfaces, which give a surprisingly large contribution to the observed strains, stresses, and swelling. It is known that in the treatment of elastic fields of point defects, strains and stresses associated with boundary conditions have the same magnitude as strains and stresses (22)–(25) generated by the defects themselves, see for example [25, 26]. Below we show that, for example, in tungsten only 59% of the observed swelling comes directly from the sources in the integral equations representing relaxation volumes of defects. The remaining 41% of the swelling effect comes from strains and stresses associated with boundary conditions.

In the next section we describe how function $\omega_{\rm rel}({\bf r})$ can be computed or estimated using ab initio methods and molecular dynamics simulations, and generic assumptions about the spatial distribution of radiation damage. We then investigate the case of a spherical shell containing an isotropic neutron source, representing for example a vacuum vessel of a fusion power plant or the casing of a fissile fuel element, and solve the above equations for the stress and strain analytically. We then develop a finite element implementation of the approach described above, to enable computing strains and stresses in components of arbitrary shape and size exposed to irradiation.

In this section we summarize results of ab initio density functional theory calculations of relaxation volumes of self-interstitial atoms (SIA) and vacancy defects in several bcc transition metals and in gold. The availability of such data is still relatively limited as the majority of density functional calculations of defects performed so far focused on the accurate evaluation of energies of defects [42–45], rather than on the evaluation of elastic properties of defects [7, 8, 27, 29, 46]. The relaxation volume of a Frenkel pair can be approximated by the sum of relaxation volumes of an SIA and a vacancy. In non-magnetic bcc transition metals, the $\langle 1\, 1\, 1 \rangle$ dumbbell is the most stable SIA defect configuration [47]. We have evaluated the relaxation volumes of SIA and vacancy defects in vanadium, niobium, molybdenum, tantalum and tungsten using equation (5). DFT calculations were performed using the approach described in [7]. Relaxation volumes of vacancies were computed using simulation boxes containing $3\times3\times3$ bcc unit cells, using a $5\times5\times5$ k-point mesh and the GGA-PBE functional.

Results of ab initio calculations of relaxation volumes of defects are summarised in table 1. The relaxation volume of a vacancy is negative whereas the relaxation volume of a SIA defect is positive, and in most cases larger than the volume of an atom. The data given in the table show that the relaxation volume of a Frenkel pair is close to one atomic volume. Relaxation volumes of point defects in fcc gold (Au) were evaluated in [41]. In gold, the most stable configuration of a SIA defect is a $\langle 1\, 0\, 0 \rangle$ dumbbell. The relaxation volume of a Frenkel pair in gold is close to 1.64 atomic volume, larger than in bcc transition metals.

Table 1. Relaxation volume of $\langle 1\, 1\, 1 \rangle$ self-interstitial atom (SIA) dumbbell defects and vacancies in non-magnetic bcc transition metals, and of a $\langle 1\, 0\, 0 \rangle$ self-interstitial dumbbell defect and a vacancy in fcc gold. The relaxation volume of a Frenkel pair is approximated by the sum of relaxation volumes of a SIA defect and a vacancy. Data are given in atomic volume units; defect relaxation volumes in gold are taken from [41].

	$\Omega^{\rm SIA}_{\rm rel}$	$\Omega^{\rm V}_{\rm rel}$	$\Omega^{\rm FP}_{\rm rel}$
V	1.472	−0.493	0.979
Nb	1.554	−0.405	1.149
Mo	1.538	−0.353	1.185
Ta	1.524	−0.410	1.114
W	1.712	−0.345	1.367
Au	2.02	−0.38	1.64

To compare defect dipole tensors computed using semi-empirical many-body interatomic potentials with ab initio results derived from density functional calculations, we evaluated dipole tensors of point defects in tungsten.

All the ab initio calculations were performed using Vienna Ab initio Simulation Package (VASP) [48–51]. We used the PBE [52, 53] and AM05 [54–56] exchange-correlation functionals. The plane wave energy cutoff was 450 eV. We used a supercell containing $4\times4\times4$ bcc unit cells, and a $5\times5\times5$ k-points mesh. First, we created perfect lattice cells containing 128 atoms and fully relaxed them, to find the equilibrium lattice parameter. Then, we created cells containing point defects. Ionic positions were relaxed, but the cell size and shape remained the same as in the perfect lattice case. Elastic dipole tensors were computed from macro-stresses using equation (5). The computed values are given in tables 2 and 3.

We have also evaluated dipole tensors of point defects using molecular statics and semi-empirical interatomic potentials. We used the Marinica (EAM4) [57] and DND potentials [47]. Similarly to ab initio calculations, we first fully relaxed a simulation box containing $80\times80\times80$ perfect bcc unit cells. Then, we inserted or removed atoms in the cell, creating point defects. This was followed by the relaxation of atomic configurations, where again we did not change the cell size and shape. The dipole tensors were computed using equation (5). The computed values are given in tables 4 and 5. Comparison of elements of dipole tensors computed using density functional theory and semi-empirical potentials in tables 2–5, show that whereas the results agree qualitatively, there are significant differences between the absolute values. This is perhaps not surprising since the main criterion used in fitting the semi-empirical potentials to ab initio data has been the agreement between the energies of atomic configurations. Elastic properties of defects have so far received relatively little attention although, as we show in this study, they play a pivotal role in determining the magnitude of strains and stresses in reactor components under irradiation.

To evaluate the density of relaxation volumes of radiation defects (21), it is not enough to know the number of Frenkel pairs per unit volume. It is only in the limit where self-interstitial atom defects and vacancies are distributed randomly as an ideal gas of defects that the density of relaxation volumes can be estimated as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \omega_{\rm rel}({\bf r})=\Omega_{\rm rel}^{\rm SIA} n_{\rm SIA}({\bf r})+\Omega_{\rm rel}^{\rm V} n_{\rm V}({\bf r}),\label{ideal_gas} \nonumber \end{align} \tag{ 31 }$$

where $n_{\rm SIA}({\bf r})$ and $n_{\rm V}({\bf r})$ are the volume densities of self-interstitial atom and vacancy defects, and where the relaxation volume of a SIA defect $\Omega_{\rm rel}^{\rm SIA}$ is positive and the relaxation volume of a vacancy $\Omega_{\rm rel}^{\rm V}$ is negative.

If defects form clusters, for example dislocation loops or voids, the relaxation volume of a cluster of defects is different from the sum of relaxation volumes of defects forming the cluster. For example, the volume of a dislocation loop in the macroscopic limit equals the scalar product of its Burgers vector and its vector area $({\bf b}\cdot {\bf A})$, and can be positive or negative, depending on the interstitial or vacancy character of the loop [7, 31]. The relaxation volume of a void, which is a large cluster of vacancies, vanishes in the macroscopic limit [39]. The relaxation volume of a vacancy-helium cluster depends on the relative number of vacancies and helium atoms in a cluster [8].

In the next section we evaluate relaxation volumes of clusters of defects formed in collision cascades. We show that the total relaxation volume of cascade debris is still proportional to the number of Frenkel pairs produced in a cascade event, however it is not given by equation (31) above.

5.1. The dipole tensor for a large, complex defect

Equation (2) gives the displacement field at a position ${\bf r}$, which is sufficiently far away from a defect so it can be treated as a point source of elastic deformation. If there are m such point sources, the displacement field at ${\bf r}$ is simply the sum over their contributions. Contributing point sources could be a substitutional atom whose local environment induces a stress upon it, but equally could be a point defect, or a cluster of defects, viewed from a distance. The dipole tensor of a large, complex defect therefore has the same simple physical interpretation as that of a point defect or a small defect cluster—it can be approximated by a point source of stress which produces the same displacement field (to leading order) as the sum of its contributing parts.

We can say that a complex defect is a group of m contributors, where the $a=\{1, 2, \ldots, m\}$th contributor has dipole tensor $P^{a}_{kl}$ and position ${\bf R}^{a}$. We would like to represent this with a single dipole tensor P_kl at position ${\bf R}$. The displacement at ${\bf r}$ due to the group of m contributors is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{displacements_complex} u_i({\bf r}) &= - \sum_{a=1}^{m}\, P^{a}_{kl} \frac{\partial}{\partial x_l} G_{ik}({\bf r} - {\bf R}^{a}) \nonumber \\ &= - \sum_{a=1}^{m}\, P^{a}_{kl} \frac{\partial}{\partial x_l} G_{ik}({\bf r - {\bf R}}) \nonumber \\ & + \sum_{a=1}^{m}\, \left(R^{a}_j - R_j\right) \, P^{a}_{kl} \, \frac{\partial^2 }{\partial x_l \partial x_j } G_{ik}({\bf r - {\bf R}}) \nonumber \\ & + \mathcal{O}\left( \frac{\left({\bf R}^{a} - {\bf R}\right)^2}{R^5} \right). \nonumber \end{align} \tag{ 32 }$$

We recognize the first term as the displacement field due to a single defect of strength $P_{kl} = \sum_{a} P^{a}_{kl}$ at ${\bf R}$, and the second term as the first-order correction due to the spatial extent of the defect. The spherical average of the second derivative of G_ik is zero, so to minimise the second term we choose

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{centreOfPosition1} {\bf R} &= {\rm argmin} \left\| \sum_{a=1}^{m}\, \left( \bf R^{a} - \bf R\right) \, P^{a}_{kl} \right\| \nonumber \\ &\simeq \frac{\sum_{a=1}^{m}\,\left\| {\bf P}^{a} \right\| {\bf R}^{a} }{\sum_{a=1}^{m}\,\left\| {\bf P}^{a} \right\| } \nonumber \end{align} \tag{ 33 }$$

where $\left\| {\bf P} \right\| = \sqrt{{\rm Tr}( {\bf P}^2) }$ is the Frobenius norm, used here as a measure of the strength of the ath contributor. In this way we can define a single dipole tensor and single position for a complex defect.

5.2. Dipole tensor for defects generated by a collision cascade

In this section we compute the dipole tensors for collision cascades in bulk tungsten simulated with classical molecular dynamics. It is not our intention here to provide an exhaustive database of values of dipole tensors, only to describe the relative ease of such a calculation. For an empirical potential we can define the energy as a sum over all atoms a, and so [58]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{kl} \equiv \sum_{a} P^{a}_{kl} = -\frac{1}{2} \sum_{a,b} f_{ab,k} r_{ab,l}, \label{dipoleTensorPerAtom} \nonumber \end{align} \tag{ 34 }$$

where f_ab,k is the kth Cartesian component of the force acting on atom a due to atom b and r_ab,l is the lth Cartesian component vector separation between them. The second derivative of energy with respect to strain gives the (fourth-rank) elastic constant tensor, needed for evaluating the relaxation volume given by equation (7).

Collision cascades were simulated with the classical molecular dynamics code PARCAS. A simulation supercell of 6.8 million atoms at 0 K was established with periodic stress-free boundary conditions, then a cascade was initiated by giving a single atom a large kinetic energy (100 keV–200 keV) in a random direction. The DND interatomic potential for tungsten by Derlet et al [47], stiffened at short range by Björkas et al [59], was used for the simulations, and a non-local friction force was applied to atoms with a kinetic energy above 10 eV to account for energy loss due to electronic stopping [60]. A Berendsen thermostat [61] set to 0 K was applied to the atoms in a 1.5 unit cell thick region along all periodic boundaries. Cascades were followed for 40 ps, at which point further defect evolution is thermal rather than ballistic. Details of the simulation method can be found in [20]. Similar cascade simulations have been shown to reproduce the observed experimental cluster size-frequency distribution [21, 23] and spatial extent [24]. A typical relaxed cascade configuration is shown in figure 2. To compute the dipole tensor for the atomically relaxed configuration, the atoms in the supercell could in principle be directly relaxed with conjugate gradients or a similar method, but in practice we have found that the dipole tensor converges rather slowly in a large box, and a standard convergence criterion based on force-per-atom or total energy difference alone may return without a well-converged dipole tensor. Instead we take the calculation as a three-stage procedure: first the atoms deviating more than $\frac{1}{4}$ lattice parameter from crystal positions were identified, and a buffer of four unit cells in all directions was added. A smaller rectangular cell containing these atoms was cropped from the cascade and fully relaxed. At this point the dipole tensor is converged to within 10%. The smaller cell was then re-embedded into a 4 million atom perfect crystal supercell ($128\times128\times128$ unit cells), and relaxed again. The dipole tensor is computed after the second relaxation. Finally the dipole tensor computed for the perfect crystal lattice at the original supercell size was computed and subtracted—this may not be exactly zero if the lattice parameter is only specified to a fixed precision, and needs removing as it scales with system size. We find convergence in the dipole tensor to better than 1% precision after a few hours cpu time. We expect this precision to be significantly better than the accuracy of the potentials or the description of the cascade structure.

Zoom In Zoom Out Reset image size

Table 6. Results for the dipole tensor for three representative displacement cascades simulated using molecular dynamics. The energy of the initial primary knock-on atom (PKA) is given in kiloelectron-Volts (keV), the number of Frenkel pairs was computed by the Wigner–Seitz cell analysis of the final cascade configuration. The partial relaxation volumes were computed using equation (8). The total relaxation volume of a cascade is also reported as a multiple of the volume per atom, $\Omega_0$. Elements of dipole tensors and relaxation volumes are converged to better than 1%.

PKA energy (keV)	100	150	200
Number of Frenkel pairs	199	174	198
Elements of the dipole tensor (eV)
P11	2340	3080	3640
P22	1480	2860	3410
P33	2030	2450	3140
P23	−907	62	−1420
P31	−1580	−335	351
P12	−367	1700	−714
Partial relaxation volumes $\Omega ^{(s)}$ ($\mathring{\rm A}^3$)
$\Omega _1$	−1570	−1150	−908
$\Omega _2$	505	331	2220
$\Omega _3$	2040	2220	2240
Relaxation volumes $\Omega _{\rm rel}/\Omega _0$	61.4	88.2	107.0

Dipole tensors were computed for isolated defects and relaxed cascade configurations using the stress method, equation (5). Results for idealised isolated defects are given in section 4, from which we conclude that empirical potentials give reasonable results. We note that the relaxation volume for a single tungsten SIA is around 1.7 $\Omega_0$, but for a large loop containing N interstitials the relaxation volume drops to $N \Omega_0$, in agreement with elasticity analysis [7, 31]. Defect dipole tensors for the full cascades are given in table 6. We conclude that the relaxation volume of a cascade scales with the number of Frenkel pairs, with the DND potential giving

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Omega_{\rm rel}(E) \simeq 0.5 N_{\rm FP}(E) \, \Omega_0. \label{relaxation_volume_cascade} \nonumber \end{align} \tag{ 35 }$$

The elastic relaxation volume density $\omega_{\rm rel}({\bf r})$ of cascades resulting from atomic recoils of energy E can be computed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \omega_{\rm rel}({\bf r})=\int n({\bf r},E)\Omega_{\rm rel}(E){\rm d}E, \nonumber \end{align} \tag{ 36 }$$

where $n({\bf r}, E)$ is the spatial density of cascades with energy E and $\Omega_{\rm rel}(E)$ is the relaxation volume of defects produced in a cascade initiated by a recoil atom of energy E. In practice, this function can be computed using neutron transport calculations for a realistic reactor geometry [63, 64] followed by the evaluation of local spectra of primary knock-on atoms [65] and the treatment of the subsequent evolution of microstructure, including modelling dislocation climb [66], self-climb [67], and the formation of a network of dislocations, voids and gas bubbles.

Above, we showed how to compute the relaxation volume density (21) using ab initio calculations and molecular dynamics simulations of high energy collision cascades. In this and the next sections we show how to solve the elasticity equations and compute strains, stresses and swelling of components assuming that function $\omega_{\rm rel}({\bf r})$ is known. In this section, we show how to solve equations (23)–(30) analytically in the limit where the problem has spherical symmetry, and in the next section we describe a general numerical finite element approach suitable for the treatment of arbitrary configurations.

Consider a component, for example the vacuum vessel of a fusion power plant or the cladding of a fuel cell containing fissile nuclear fuel, which we assume has the form of a spherical shell sketched in figure 3.

Zoom In Zoom Out Reset image size

If the source of neutrons inside the shell is isotropic, then the distribution of defect relaxation volumes in the material depends only on the distance r to the centre of the shell and is independent of the polar and azimuthal angles θ and ϕ of the spherical system of coordinates, the origin of which is at the centre of the shell.

The solution that we are interested in is defined on the interval $R_1\leqslant r \leqslant R_2$, where R₁ is the inner radius of the shell and R₂ is its outer radius. The field of displacements is radially-symmetric and the vector of displacements can be written as ${\bf u}({\bf r})=u_r(r){\bf n}$, where ${\bf n}={\bf r}/r$. Differentiating equation (30), we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\partial \over \partial x_i} u_i({\bf r})&=-{1\over 12\pi}\left({1+\nu \over 1-\nu}\right){\partial^2 \over \partial x_i^2}\int {\omega _{\rm rel}({\bf R}) \over |{\bf r}-{\bf R}|} {\rm d}^3R\nonumber \\ &={1\over 12\pi}\left({1+\nu \over 1-\nu}\right)\int \omega _{\rm rel}({\bf R}) [4\pi \delta({\bf r}-{\bf R})] {\rm d}^3R\nonumber \\ &= {1\over 3}\left({1+\nu \over 1-\nu}\right)\omega_{\rm rel}({\bf r}) .\label{displacements_divergence} \nonumber \end{align} \tag{ 37 }$$

The above equation has a simple meaning, namely that it is the density of relaxation volumes of defects that causes the material to expand or contract as a result of accumulation of defects in the material. Note the remarkable prefactor $(1/3)(1+\nu)/(1-\nu)$ in the above equation, the numerical value of which is close to 0.6 for tungsten. This prefactor shows that the direct deformation of the lattice caused by the accumulation of defects is only partially responsible for the observed dimensional changes occurring as a result of irradiation. The other, similar in magnitude, contribution comes from the boundary conditions at surfaces, as we prove below.

Since the field of displacements is radially symmetric, we apply the divergence theorem to equation (37) and write

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 4\pi r^2u_r(r)={4\pi\over 3}\left({1+\nu \over 1-\nu}\right)\int _{R_1}^r R^2\omega _{\rm rel}(R){\rm d}R,\label{u_r} \nonumber \end{align} \tag{ 38 }$$

for $R_1\leqslant r\leqslant R_2$.

Noting that in the absence of radiation defects the divergence of ${\bf u}({\bf r})$ is a harmonic function of coordinates [40], we write the field of displacements as a sum of the partial solution of heterogeneous equation (37) and a general solution of the corresponding homogeneous equation, namely

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u_r(r)=ar+{b\over r^2}+{1\over 3r^2}\left({1+\nu \over 1-\nu}\right)\int _{R_1}^r R^2\omega _{\rm rel}(R){\rm d}R.\label{u_r_full} \nonumber \end{align} \tag{ 39 }$$

Here a and b are constants that need to be determined from boundary conditions at r  =  R₁ and r  =  R₂. Assuming that pressure at R₁ and R₂ is negligible in comparison with the stresses developing in the material due to the accumulation of defects, we adopt the traction free boundary conditions [25] $\sigma _{ij}(r)n_j=0$ at r  =  R₁ and r  =  R₂.

Strains can be found by differentiating (39). They have the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon _{ij}(r) =&~ a\delta_{ij} + {b\over r^3}\left(\delta _{ij}-3n_in_j\right)\nonumber \\ &+ {1\over 3r^3}\left({1+\nu \over 1-\nu}\right) \left(\delta _{ij}-3n_in_j\right)\int _{R_1}^r R^2\omega _{\rm rel}(R){\rm d}R\nonumber \\ &+{1\over 3}\left({1+\nu \over 1-\nu}\right) n_i n_j \omega_{\rm rel}({r}). \label{strain_full} \nonumber \end{align} \tag{ 40 }$$

To find stresses, we need to multiply the above expression for the strain tensor by the four-index tensor of elastic constants (19). The resulting expression has the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{ij}(r)=&~2\mu a\left({1+\nu \over 1-2\nu}\right)\delta_{ij}+{2\mu b\over r^3}(\delta_{ij}-3n_in_j)\nonumber \\ &+{2\mu \over 3r^3}\left({1+\nu \over 1-\nu}\right)\left(\delta _{ij}-3n_in_j\right)\int _{R_1}^r R^2\omega _{\rm rel}(R){\rm d}R\nonumber \\ &+{\mu\over 3}\left({1+\nu \over 1-\nu}\right) \omega_{\rm rel}({r})\left( {2\nu \over 1-2\nu}\delta_{ij}+2 n_in_j\right). \label{stress_full} \nonumber \end{align} \tag{ 41 }$$

Applying the traction-free boundary conditions at r  =  R₁ and r  =  R₂ to (41), we find parameters a and b, namely

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle a={1\over 6\pi}\left({1-2\nu \over 1-\nu}\right){\Omega_{\rm tot}\over R_2^3-R_1^3},\label{a} \nonumber \end{align} \tag{ 42 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle b={1\over 12\pi}\left({1+\nu \over 1-\nu}\right){R_1^3\over R_2^3-R_1^3}\Omega_{\rm tot},\label{b} \nonumber \end{align} \tag{ 43 }$$

where

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \Omega _{\rm tot}=4\pi \int _{R_1}^{R_2}R^2\omega_{\rm rel}(R){\rm d}R \nonumber \end{align*}$$

is the total relaxation volume of all the defects in the material of the shell.

The total macroscopic change of volume resulting from the accumulation of defects in the materials of the shell is given by the integral of the trace of the full strain tensor (40) over the volume of the component

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta V&=\int \epsilon _{ii}({\bf r}){\rm d}V=\int \left[{1\over 3}\left({1+\nu \over 1-\nu}\right) \omega_{\rm rel}({\bf r})+3a\right]{\rm d}V\nonumber \\ &=\int _{R_1}^{R_2}\left[{1\over 3}\left({1+\nu \over 1-\nu}\right) \omega_{\rm rel}(r)+3a\right]4\pi r^2{\rm d}r, \label{swelling} \nonumber \end{align} \tag{ 44 }$$

where we noted that $\delta_{ii}=3$ and $n_in_i=1$. The integration of the first term in square brackets over the volume of the shell gives

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \Delta V_1={1\over 3}\left({1+\nu \over 1-\nu}\right)\Omega_{\rm tot}, \nonumber \end{align*}$$

whereas the second term in (44), proportional to a and arising from the boundary conditions, contributes

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \Delta V_2={2\over 3}\left({1-2\nu \over 1-\nu}\right) \Omega_{\rm tot} \nonumber \end{align*}$$

to the total macroscopic swelling of the shell. The sum of these two terms equals $\Omega_{\rm tot}$, confirming that swelling is as much an effect of direct expansion of the lattice due to the accumulation of defects in the material, as it is an effect associated with boundary conditions and arising from long-range stresses and strains of defects interacting with the boundaries of the component.

Substituting (42) and (43) into (39), we find radial displacements of the inner and outer surfaces of the shell

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle u_r(R_1)&={1\over 4\pi}{R_1\over R_2^3-R_1^3}\Omega_{\rm tot}\nonumber \\ u_r(R_2)&={1\over 4\pi}{R_2\over R_2^3-R_1^3}\Omega_{\rm tot}.\label{inner_and_outer_surfaces} \nonumber \end{align} \tag{ 45 }$$

These displacements satisfy the condition

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle 4\pi R^2_2u_r(R_2)-4\pi R_1^2u_r(R_1)=\Omega_{\rm tot},\label{swelling_surfaces} \nonumber \end{align} \tag{ 46 }$$

which provides an alternative way of evaluating the total volumetric swelling of the shell. A remarkable property of equation (45) is that they show that both the inner and outer surfaces of an irradiated spherical shell relax outwards. The magnitude of surface displacements can be substantial. For example, the surface of a vacuum vessel with the inner radius of R₁  =  5 metres, after exposure to neutron irradiation, producing constant homogeneous defect relaxation volume density $\omega_{\rm rel}=1$% in the material, moves outwards by approximately $u_r(R_1)=\omega_{\rm rel}R_1/3\approx 1.6$ cm.

We can now find elements of the stress tensor developing in the spherical shell as a result of its exposure to irradiation. The radial diagonal element of the stress tensor is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{rr}(r)=&~{\mu \over 3\pi}\left({1+\nu \over 1-\nu}\right)\left[{\Omega_{\rm tot}\over R_2^3-R_1^3 }\left(1-{R_1^3\over r^3}\right)\right]\nonumber \\ &-{\mu \over 3\pi}\left({1+\nu \over 1-\nu}\right)\left[{4\pi\over r^3}\int _{R_1}^r R^2\omega_{\rm rel}(R){\rm d}R\right]\nonumber \\ &+{2\mu\over 3}\left({1+\nu \over 1-2\nu}\right) \omega_{\rm rel}(r).\label{stress_rr} \nonumber \end{align} \tag{ 47 }$$

The circumferential (hoop) components $\sigma_{\theta\theta}(r)$ and $\sigma_{\phi\phi}(r)$ of the stress tensors are

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{\theta\theta}(r)=&~{\mu \over 3\pi}\left({1+\nu \over 1-\nu}\right)\left[{\Omega_{\rm tot}\over R_2^3-R_1^3 }\left(1+{R_1^3\over 2r^3}\right)\right]\nonumber \\ &+{\mu \over 6\pi}\left({1+\nu \over 1-\nu}\right)\left[{4\pi\over r^3}\int _{R_1}^r R^2\omega_{\rm rel}(R){\rm d}R\right]\nonumber \\ &+{\mu\over 3}\left({1+\nu \over 1-\nu}\right) \left({2\nu \over 1-2\nu}\right)\omega_{\rm rel}(r).\label{stress_ff} \nonumber \end{align} \tag{ 48 }$$

Due to the spherical symmetry of the problem, we have $\sigma_{\theta\theta}(r)=\sigma_{\phi\phi}(r)$. The above formulae are valid for any distribution of defect relaxation volume density $\omega _{\rm rel}(r)$ in the spherical shell, which is a function of radial variable r.

As the density of defects vanishes in the immediate vicinity of free surfaces, $\omega_{\rm rel}(R_1)=\omega_{\rm rel}(R_2)=0$. The radial stress (47) also vanishes at surfaces at r  =  R₁ and r  =  R₂, but the hoop stress (48) remains finite everywhere on the interval $R_1\leqslant r \leqslant R_2$. The high radial and hoop stresses developing as a result of accumulation of defects in the material are responsible for the loss of structural integrity of the component in the limit where the relaxation volume density of defects exceeds a certain critical level.

Applications of the above equations to the evaluation of stresses developing in a steel vacuum vessel with inner radius R₁  =  3 m and outer radius R₂  =  3.5 m as a result of irradiation are illustrated in figure 4. The form of function $\omega_{\rm rel}(r)$ used as input reflects the fact that neutron flux from an isotropic source varies as a function of r as $\omega_{\rm rel}(r)\sim r^{-2}$ and also that neutrons are absorbed in steel over a characteristic distance of a fraction of a metre.

Zoom In Zoom Out Reset image size

Figure 4 shows that particularly high stresses develop close to the internal surface of the steel shell even if the amount of damage accumulated in the material is relatively small. The magnitude of radiation-induced stresses can be estimated by evaluating the product of the shear modulus of the material and function $\omega_{\rm rel}({\bf r})$, suitably averaged over the volume of the component.

In the preceding section, to find analytical solutions for the displacements, strains and stresses in a hollow spherical shell exposed to irradiation, we used equation (30) as a starting point for the analysis. In this section, we develop a general finite element implementation of the method, using which one can evaluate stresses, strains and displacements in any reactor geometry. In the finite element method (FEM), it is more convenient to use the equation for body forces generated by defects (29) rather than the formula for displacements (30). The two approaches are entirely equivalent, as confirmed by the comparison of numerical solutions derived using analytical and FEM implementations for the same reactor geometry and same density of relaxation volumes of defects (see figure 4).

The finite element implementation of the mathematical formalism for computing strains, stresses and swelling of irradiated reactor components developed above, is based on equation (29), which can be written in the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_i({\bf r})=-B{\partial \over \partial x_i}\omega_{\rm rel}({\bf r}),\label{body_force2} \nonumber \end{align} \tag{ 49 }$$

where $B=2\mu(1+\nu)/[3(1-2\nu)]$ is the bulk modulus of the material and $\omega_{\rm rel}({\bf r})$ is the density of relaxation volumes of defects (21).

The fundamental equation for the stress tensor, defining the condition of mechanical equilibrium, has the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\partial \over \partial x_j}\sigma _{ij}({\bf r})+f_i({\bf r})=0.\label{condition_of_equilibrium} \nonumber \end{align} \tag{ 50 }$$

To find the displacement field which satisfies (50) using the FEM, this expression for mechanical equilibrium is recast as a virtual work expression [69]; for completeness this is outlined here. Multiplying (50) by an arbitrary displacement field $\delta u_i$ which satisfies $\delta u_i = 0$ on the displacement boundary S_u, and integrating over the body $V$ gives:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:weak} \int_V \left({\partial\sigma _{ij} \over \partial x_j}+f_i\right)\delta u_i {\rm d}^3R=0. \nonumber \end{align} \tag{ 51 }$$

Consider the first term

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_V {\partial \sigma_{ij} \over \partial x_j} \delta u_i {\rm d}^3R &= \int_V {\partial \over \partial x_j}\left(\sigma_{ij} \delta u_i \right) {\rm d}^3R \nonumber \\ &- \int_V \sigma_{ij} {\partial \over \partial x_j} \delta u_i {\rm d}^3R. \nonumber \end{align} \tag{ 52 }$$

Using the divergence theorem, the first term on the right hand side is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_V {\partial \over \partial x_j}\left(\sigma_{ij} \delta u_i \right) {\rm d}^3R &= \int_{S} \sigma_{ij}\delta u_i n_j {\rm d}^2R \nonumber \\ &= \int_{S_t} \delta u_i \bar{t}_i {\rm d}^2R, \nonumber \end{align} \tag{ 53 }$$

where $\bar{t}_i = \sigma_{ij}n_j$ are the specified surface tractions on S_t and by definition the virtual displacement $\delta u_i = 0$ on S_u and so the surface integral vanishes there. Therefore, after multiplying by  −1, (51) becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_V \sigma _{ij} {\partial\over \partial x_j}\delta u_i {\rm d}^3R - \int_V f_i\delta u_i {\rm d}^3R - \int_{S_t}\bar{t}_i\delta u_i {\rm d}^2R = 0. \nonumber \end{align} \tag{ 54 }$$

Expressing the stress tensor in terms of displacements,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sigma_{ij} = {1 \over 2}C_{ijkl}\left({\partial u_k\over \partial x_l}+{\partial u_l\over \partial x_k} \right) = C_{ijkl}{\partial u_k\over \partial x_l}, \nonumber \end{align} \tag{ 55 }$$

due to the symmetry of C_ijkl. Therefore

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:virtual_work} &\int_V C_{jikl}{\partial u_k\over \partial x_l} {\partial \delta u_i\over \partial x_j} {\rm d}^3R - \int_V f_i\delta u_i {\rm d}^3R \nonumber \\ &- \int_{S_t}\bar{t}_i\delta u_i {\rm d}^2R = 0. \nonumber \end{align} \tag{ 56 }$$

This is the expression for virtual work in the presence of a body force f_i. This can be discretised and solved at the FE nodes located at $r_i^a$ where 1  <  a  <  m, for the unknown displacement at each node $U_i^a$. The virtual displacement and its derivative are interpolated using the shape functions

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:deltau} \delta u_i({\bf r}) = \sum_{a=1}^m N^a({\bf r})\delta U_i^a, \nonumber \end{align} \tag{ 57 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\partial \delta u_i({\bf r}) \over \partial x_j} = \sum_{a=1}^m {\partial N^a({\bf r}) \over \partial x_j}\delta U_i^a. \nonumber \end{align} \tag{ 58 }$$

The displacement derivative is interpolated in a similar manner, namely

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\partial u_k({\bf r}) \over \partial x_l} = \sum_{b=1}^m {\partial N^b({\bf r}) \over \partial x_l}U_k^b \label{eq:dudx}, \nonumber \end{align} \tag{ 59 }$$

where $N^a({\bf r})$ is the nodal shape function of node a which satisfies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle N^a({\bf r}^b) = \delta_{ab}. \nonumber \end{align} \tag{ 60 }$$

Substituting (57)–(59) into (56) gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &\left[\int_V C_{ijkl}{\partial N^b({\bf r}) \over \partial x_l}U_k^b {\partial N^a({\bf r}) \over \partial x_j} {\rm d}^3R \right. \nonumber \\ &- \int_V f_i({\bf r}) N^a({\bf r}) {\rm d}^3R - \left.\int_{S_t} \bar{t}_i({\bf r}) N^a({\bf r}) {\rm d}^2R \right]\delta U_i^a = 0, \nonumber \end{align} \tag{ 61 }$$

and as the virtual displacement $\delta U_i^a$ is arbitrary, the term in square brackets must be zero. In matrix form, the following system of equations is obtained

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}ll@{}} K_{ik}^{ab} U_k^b = F_i^a & {\rm if }~~r_k^b~~{\rm not~~on }~~S_u \nonumber \\ U_k^b = \bar{u}_k({\bf r}^b) & {\rm if }~~r_k^b~~{\rm on }~~S_u, \end{array} \label{eq:FEM} \nonumber \end{align} \tag{ 62 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle K_{ik}^{ab} &= \int_V C_{ijkl}\frac{\partial N^a({\bf r})}{\partial x_j}\frac{\partial N^b({\bf r})}{\partial x_l} {\rm d}^3R\nonumber \\ F_i^a &= \int_V f_i({\bf r}) N^a({\bf r}) {\rm d}^3R + \int_{S_t} \bar{t}_i({\bf r}) N^a({\bf r}) {\rm d}^2R. \nonumber \label{eq:F} \nonumber \end{align} \tag{ 63 }$$

$K_{ik}^{ab}$ is the global stiffness matrix. $F_i^a$ is the global force vector which includes a body force contribution at every node not on S_u, and surface traction contribution $\bar{t}_i$ for every node on S_t. Substituting the body force $f_i({\bf r})$ from (49) into (63) gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle F_i^a = -B\int_V {\partial\omega_{\rm rel}({\bf r}) \over \partial x_i}N^a({\bf r}) {\rm d}^3R + \int_{S_t} \bar{t}_i({\bf r}) N^a({\bf r}) {\rm d}^2R. \nonumber \end{align} \tag{ 64 }$$

The advantage of using the body force is that it has a remarkably simple form (49). Furthermore any desired traction and displacement boundary conditions on the surface can be applied directly without modification. If $\omega_{\rm rel}(r)$ is an analytic function which can be differentiated then the global force vector $F_i^a$ can be specified exactly, and the nodal displacements obtained by solving [69] (62). If only numerical values $\omega_{\rm rel}^a$ are known at the nodal positions ${\bf r}^a$, then the gradient of the density of defect relaxation volumes required to obtain the body force can instead be evaluated using the finite element shape functions, in the same manner as for the displacements,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\partial\over \partial x_i}\omega_{\rm rel}({\bf r}) = \sum_{a=1}^n {\partial N^a({\bf r})\over\partial x_i} \omega_{\rm rel}^a. \nonumber \end{align} \tag{ 65 }$$

A body force can easily be implemented in a commercial finite element code such as Abaqus [70] through the DLOAD user subroutine to specify $f_i({\bf r})$ at the integration points. The finite element implementation was first validated by simulating the spherical shell solved analytically in section 6. Due to the symmetry of the problem only one quarter of the shell was simulated with symmetry boundary conditions applied. An elastic isotropic material law was used with $\nu=0.29$ and $\mu = 80$ GPa. A convergence study was performed, and $2.55\times 10^6$ antisymmetric elements (CAX4) of 1 mm in size were found to produce excellent agreement with the analytic solution, as shown in figure 4.

Since the FEM implementation is not constrained by any symmetry considerations and can be applied to compute stresses and strains for any distribution of radiation defects, here we also performed FEM simulations for a case where $\omega({\bf r})$ is no longer a function of only radial variable r but also depends on the polar angle θ. We consider the density of relaxation volumes of the form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \omega_{\rm rel}({\bf r}) = \omega_{\rm rel}(r,\theta)=\left\{\begin{array}{@{}lc@{}} 0, & \quad {\rm if }~~0\leqslant\theta \leqslant \pi/9 \nonumber \\ \omega_{\rm rel}(r), & \quad {\rm if }~~\pi/9< \theta \leqslant \pi. \end{array} \right. \nonumber \end{align} \tag{ 66 }$$

This density of relaxation volumes of defects represents the case where the material in the shell is shielded from neutrons in a segment of the solid angle within $20^\circ$ of the upper pole of the spherical shell. As the body force is no longer spherically symmetric and explicitly depends on θ, $\sigma_{\theta\theta}$ now differs from $\sigma_{\phi\phi}$ as shown in figure 5. Also, the stress field now has a non-vanishing $\sigma_{r\theta}(r, \theta)$ component, the magnitude of which is the largest at the boundary of the shielded region corresponding to $\theta =\pi/9$. The problem still has the axial symmetry as the density of relaxation volumes of defects is independent of angle ϕ.

Zoom In Zoom Out Reset image size

Half of the spherical shell was simulated at sufficiently high resolution, requiring $5.10\times10^6$ elements. Unlike in the spherically symmetric cases 1 and 2 illustrated in figure 4, the shear stress $\sigma_{r\theta}$ is now non-zero. Near the interface between the shielded and unshielded regions a stress concentration of $\sigma _{r\theta}$ is clearly visible, as shown in figure 6. Away from the the shielded region the stress is very close to the the analytical solution and the symmetric FEM solution; this is as expected. This case study demonstrates that the FEM implementation of the method based on equations (49) and (50) is completely general and can be applied to any reactor geometry and boundary conditions, provided a sufficiently fine FE mesh is used.

Zoom In Zoom Out Reset image size

The treatment developed in this paper shows how to compute stresses, strains and swelling in components of a nuclear fusion or fission reactor resulting from the accumulation of defects in materials exposed to irradiation. By deriving the fundamental macroscopic equations from the atomic scale, we show that populations of radiation defects accumulating in the materials produce internal body forces (29) and (49) that result in the build-up of internal stresses and cause macroscopic deformation of components. The central notion of the treatment is the spatially varying density of relaxation volumes of defects (21), which acts as a source of stresses and strains, and the gradient of which generates the spatially distributed body force (29). We illustrate applications of the method using exact analytical solutions and numerical finite element implementations, which suggest that the method is able to provide a reasonably accurate first-principles based foundation for the in silico engineering assessment of operational performance of a fusion power plant and its design.

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 under grant agreement No 633053 and from the RCUK Energy Programme [grant number EP/P012450/1]. The views and opinions expressed herein do not necessarily reflect those of the European Commission. We gratefully acknowledge stimulating discussions with A.P. Sutton, C.H. Woo, F. Hofmann, M. Rieth, M. O'Brien and J.W. Connor. SLD is grateful to I. Lindsey and M.M.E. Coolsen who enabled this research to be carried out. ET acknowledges support from EPSRC through fellowship grant, EP/N007239/1. AES acknowledges support from the Academy of Finland through project No. 311472.