Macroscopic elastic stress and strain produced by irradiation

The accumulation of microscopic defects resulting from the exposure of materials to neutron irradiation is one of the critical factors limiting the lifetime of components in a fission power plant [1] and, by implication, is an issue of fundamental significance in the context of a comparative assessment of designs of power-generating fusion reactors [2, 3]. The ambitious timescale for the development of fusion reactors forces one to rely on virtual engineering [4, 5], hence creating a strong drive towards the development of the so-called 'digital twins', a relatively new concept in the area of nuclear fusion [6]. It is important to include in this process the consideration of operating conditions, projecting the reactor design beyond the digital twin concept, to enable the assessment of effects of neutron irradiation. This point was highlighted by Jumel et al [7] in relation to the problem of ageing of fission reactors. The present study and our earlier work [8] are the initial steps in establishing the fundamental link between the microscopic effects of exposure of materials to irradiation and macroscopic component-scale finite element virtual reactor models.

In a fission reactor, cladding materials receive the highest exposure to neutron irradiation [9, 10], comparable to the radiation exposure expected in a fusion power plant. The accumulation of defects in a zirconium alloy, a material often used as cladding for uranium oxide fissile nuclear fuel, gives rise to anisotropic dimensional changes and irradiation-induced deformations in the limit of high irradiation dose. Since zirconium alloys have hexagonal close-packed (hcp) crystal structure, the observed dimensional changes stem primarily from the texture of the cladding material [10, 11] and the intrinsic anisotropy of the hcp crystal structure. Not only the point defects in hcp crystals are anisotropic [12], but the dynamics of agglomeration of defects into mesoscopic and macroscopic dislocation loops and dislocation networks is anisotropic as well [10, 13], resulting in significant macroscopic deformations developing in the limit of high irradiation dose [9, 10].

On the other hand, the candidate materials typically selected for the structural components of a fusion power plant all have isotropic cubic crystal structure. It is either the body-centred cubic structure characterising iron–chromium alloys, ferritic steels and tungsten, or the face-centred cubic structure of copper, copper alloys and austenitic steels [2, 14–17]. Exposure of these materials to irradiation gives rise to swelling or, in other words, to the predominantly isotropic volume expansion resulting from the accumulation of defects in the material [18–21]. In the absence of applied load or other constraints, no significant anisotropic dimensional changes are expected to occur in the crystallographically isotropic materials selected for the structural components of a fusion power plant.

A remaining source of anisotropy stems from the macroscopic spatial variation of the neutron exposure itself. Components of a fusion power plant are going to be exposed to neutron and ion irradiation generated by a localised source, the fusion plasma. The flux of high energy 14.1 MeV neutrons attenuates in the blanket surrounding the plasma on the length scale of several tens of centimetres [22, 23]. Energetic ions penetrate into the plasma-facing materials to the depth of several micrometres [24–26], whereas hydrogen and helium produced either by transmutation or by a direct exposure to the fusion plasma, can diffuse deep into the bulk of structural components, particularly at elevated temperatures [27].

Anisotropic stress in reactor components stems not only from the formation of radiation defects. Gravitational body forces act in the direction of g and the magnitude of gravitational stress is appreciable, for example the 54 divertor cassettes in ITER weighs approximately 10 tons each [28]. The magnetic field in a tokamak exerts magnetic body forces on the ferromagnetic ferritic–martensitic steels in the breeding blanket modules, and these forces are expected to be several times stronger than gravity. Thermal expansion is also expected to generate high stresses in reactor components [29]. The difference between the forces from gravity, magnetic fields and thermal expansion, and the internal forces from the defects formed following the exposure of materials to neutron irradiation is that irradiation changes the microscopic structure of a material, and hence the effects of irradiation, including radiation-induced deformations and stress, are in most cases irreversible.

The formation of radiation defects [30, 31] and defects associated with gas impurity atoms [32] give rise to the volumetric expansion of the material [33]. The magnitude of this volumetric expansion varies depending on the irradiation dose and temperature [20], whereas the presence of helium and hydrogen may have a significant additional non-linear effect on the microstructure and swelling rates [34–36].

Swelling typically increases monotonically with exposure to irradiation and exhibits a complex pattern of variation as a function of dose and temperature [20, 37]. It is natural to expect that in an operating reactor, where different parts of every component have different temperatures and are exposed to irradiation at different dose rates, swelling is spatially non-uniform. In this study, we investigate what effect this spatial non-uniformity of swelling is going to have on the stress and strain developing in reactor components during the operation of the reactor. The non-uniform irradiation, temperature and loading conditions developing in large-scale reactor components are different from the uniform or zero-stress state developing in test samples of materials employed in the context of irradiation experiments and tests [15, 38].

In an earlier study [8], we found that the density of relaxation volumes of defects represents the fundamental notion linking microscopic and macroscopic scales in a holistic simulation of radiation-induced deformations. The treatment developed in [8] was focussed solely on the evaluation of the total stress and strain fields, which include both the elastic and defects' own contributions to deformation.

However, only the elastic stress enters the engineering structural integrity criteria. Here, we develop an approach that enables computing the purely elastic stresses and deformations in an irradiated component on the macroscale, with an additional emphasis on the effect of boundary conditions.

The central concept that enables separating the purely elastic part of deformation from the total deformation, and at the same time developing a genuinely multiscale approach to the problem, is the notion of Mura's eigenstrain [39], and the fact that it can be identified uniquely with the tensorial density of relaxation volumes of defects. This statement and a proof of its validity is given in the first section of the paper.

We then explore several case studies involving representative geometries of components or samples exposed to irradiation. They include a specimen of rectangular geometry often used in the context of ion irradiation experiments. We then explore the deformations developing in a pressurised cylindrical pipe exposed to irradiation from an external source. This case describes a typical test sample used for investigating the combined effect of stress and neutron irradiation on materials. The most accessible to mathematical analysis case is a spherical shell exposed to a centrally-symmetric source of irradiation. This case admits an exact analytical solution, which agrees fully with numerical finite element model (FEM) simulations. Finally, using a purely numerical FEM approach, we evaluate elastic stress and strain developing in a cellular rectangular component exposed to irradiation. Such a cellular structure resembles a section of a tritium breeding blanket of a fusion tokamak reactor, and the simulations illustrate a complex pattern of stress and strain that is expected to develop in a geometrically complex macroscopic reactor component exposed to a spatially varying neutron flux—entirely different from the cases where irradiation is more uniform.

In this work, we adopt the same approximation as in reference [8], where the density of relaxation volume tensors of defects is taken as diagonal with respect to its Cartesian components. This spatially-varying isotropic volumetric swelling approximation is valid if the defects adopt random configurations with respect to their internal directional degrees of freedom [8]. There are cases where the applied stress, or the stress generated by the defects themselves, polarize the relaxation volume tensor density of defects [25, 40, 41]. This aspect of microstructural evolution and its effect on the elastic stress and strain fields, while being accessible to the methodology developed below, is a non-linear microstructural relaxation phenomenon that remains outside the scope of this study.

When investigating the various case studies below, we take the local relaxation volume density of defects at a given location as proportional to the radiation exposure of the material at that location. In a real operating scenario, the density of relaxation volumes at a given location is given by a complex convolution of histories of radiation exposure, temperature and stress at that location, see for example references [25, 42], which themselves depend in a self-consistent manner on the history of evolution of these quantities at other locations in a component and possibly even in the entire reactor structure. Still, irrespectively of the approximations adopted for the evaluation of the density of relaxation volumes of defects, the fundamental link between this entity and the elastic strain and stress fields generated by irradiation remains unique, and it is this link and its implications that we establish and explore in detail below.

Since the spatial scale of a macroscopic component is many orders of magnitude larger than the microscopic size of defects present in the material, the field of displacements u(x) generated by a defect situated at x' can be evaluated using the far-field approximation as

$$\begin{equation}{u}_{i}(\mathbf{x})=-{P}_{jl}\frac{\partial }{\partial {x}_{l}}{G}_{ij}(\mathbf{x}-{\mathbf{x}}^{\prime }),\end{equation} \tag{ 1 }$$

where P_jl is the elastic dipole tensor of the defect [43, 44], and G_ij (x − x') is the Green function of elasticity equations [39, 45]. A general form of equation (1), describing the field of displacements produced by an ensemble of defects distributed in the material, is

$$\begin{equation}{u}_{i}(\mathbf{x})=-\int {{\Pi}}_{jl}({\mathbf{x}}^{\prime })\frac{\partial }{\partial {x}_{l}}{G}_{ij}(\mathbf{x}-{\mathbf{x}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime },\end{equation} \tag{ 2 }$$

where ${{\Pi}}_{jl}(\mathbf{x})={\sum }_{a}{P}_{jl}^{(a)}\delta (\mathbf{x}-{\mathbf{R}}_{a})$ is the density of elastic dipole tensors of defects. This field is related to another field ω_mn (x), describing the density of relaxation volume tensors of defects, via

$$\begin{equation}{{\Pi}}_{jl}(\mathbf{x})={C}_{jlmn}{\omega }_{mn}(\mathbf{x}),\end{equation} \tag{ 3 }$$

where C_jlmn is the fourth-rank tensor of elastic constants. The density of relaxation volume tensors ω_mn (x) is defined as

$$\begin{equation}{\omega }_{mn}(\mathbf{x})=\sum\limits _{a}{{\Omega}}_{mn}^{(a)}\delta (\mathbf{x}-{\mathbf{R}}_{a}),\end{equation} \tag{ 4 }$$

where ${{\Omega}}_{mn}^{(a)}$ is the relaxation volume tensor of a defect situated at R_a . Relaxation volume tensors of defects can now be routinely computed using density functional theory [30, 31, 46, 47] or, with lesser accuracy, using molecular dynamics simulations, where interactions between the atoms are described by semi-empirical potentials [33, 48]. For a particular defect configuration, the tensorial density of relaxation volumes is given by equation (4), whereas in the general case it can be treated as an ensemble-averaged quantity [8]

$$\begin{equation}{\omega }_{mn}(\mathbf{x})=\left\langle \sum\limits _{a}{{\Omega}}_{mn}^{(a)}\delta (\mathbf{x}-{\mathbf{R}}_{a})\right\rangle ,\end{equation} \tag{ 5 }$$

where ⟨⋯⟩ refers to averaging over a statistical ensemble of representative realisations of microstructure. Function (4) is singular at locations occupied by the defects, whereas (5) is non-singular.

In terms of the density of relaxation volume tensors of defects (4) and (5), the field of displacements can be expressed as

$$\begin{equation}{u}_{i}(\mathbf{x})=-\int {C}_{jlmn}{\omega }_{mn}({\mathbf{x}}^{\prime })\frac{\partial }{\partial {x}_{l}}{G}_{ij}(\mathbf{x}-{\mathbf{x}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime }.\end{equation} \tag{ 6 }$$

Comparing this relation with equation (3.23) by Mura [39]

$$\begin{equation}{u}_{i}(\mathbf{x})=-\int {C}_{jlmn}{{\epsilon}}_{mn}^{\ast }({\mathbf{x}}^{\prime })\frac{\partial }{\partial {x}_{l}}{G}_{ij}(\mathbf{x}-{\mathbf{x}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime },\end{equation} \tag{ 7 }$$

we establish and prove the defect eigenstrain theorem, namely the statement that Mura's spatially-varying field of eigenstrains ${{\epsilon}}_{mn}^{\ast }(\mathbf{x})$ is identical to the density of relaxation volume tensors of defects, viz

$$\begin{equation}{{\epsilon}}_{mn}^{\ast }(\mathbf{x})\equiv {\omega }_{mn}(\mathbf{x}).\end{equation} \tag{ 8 }$$

The recognition of equivalence of eigenstrain ${{\epsilon}}_{mn}^{\ast }(\mathbf{x})$ to the density of relaxation volume tensors ω_mn (x) enables defining pure elastic strain and stress following the convention [39]

$$\begin{equation}{{\epsilon}}_{ij}(\mathbf{x})={{\epsilon}}_{ij}^{(\text{tot})}(\mathbf{x})-{\omega }_{ij}(\mathbf{x}),\end{equation} \tag{ 9 }$$

and

$$\begin{align}\hfill {\sigma }_{ij}(\mathbf{x})& ={C}_{ijkl}{{\epsilon}}_{kl}(\mathbf{x})={\sigma }_{ij}^{(\text{tot})}(\mathbf{x})-{{\Pi}}_{ij}(\mathbf{x})\hfill \\ \hfill & ={\sigma }_{ij}^{(\text{tot})}(\mathbf{x})-{C}_{ijkl}{\omega }_{kl}(\mathbf{x}),\hfill \end{align} \tag{ 10 }$$

where ${\sigma }_{ij}^{(\text{tot})}(\mathbf{x})={C}_{ijkl}{{\epsilon}}_{kl}^{(\text{tot})}(\mathbf{x})$. The total strain in (9) is taken as a symmetrised derivative of atomic displacements [49]

$$\begin{equation}{{\epsilon}}_{kl}^{(\text{tot})}(\mathbf{x})=\frac{1}{2}\left(\frac{\partial {u}_{k}(\mathbf{x})}{\partial {x}_{l}}+\frac{\partial {u}_{l}(\mathbf{x})}{\partial {x}_{k}}\right),\end{equation} \tag{ 11 }$$

and hence equations (9) and (10) define the purely elastic components of spatially varying strain and stress. Equations (9) and (10) show that at the macroscopic level, elastic strain and stress fields are defined taking the notion of a material containing a certain local concentration of defects, corresponding to eigenstrain ω_ij (x), as a reference non-distorted configuration.

Defects can be produced not only by irradiation but also by processing, for example for the purpose of creating a microstructurally complex material. A complex microstructure can be viewed as a highly imperfect material, characterised by a high density of regions where local atomic structure deviates significantly from ideal crystalline order. An example of a microstructurally complex material is steel, where strongly fluctuating plastic strains and deformations are introduced during manufacturing by alloying and through thermo-mechanical treatments.

In the absence of conventional macroscopic body forces, for example gravity or thermal expansion [49], the condition of mechanical equilibrium [45] ∂σ_ij (x)/∂x_j = 0, expressed in terms of elastic stress (10), has the form

$$\begin{equation}{C}_{ijkl}\frac{\partial {{\epsilon}}_{kl}^{(\text{tot})}(\mathbf{x})}{\partial {x}_{j}}-{C}_{ipmn}\frac{\partial {\omega }_{mn}(\mathbf{x})}{\partial {x}_{p}}=0.\end{equation} \tag{ 12 }$$

Interpreting this as a condition of equilibrium for the total stress, namely $\partial {\sigma }_{ij}^{(\text{tot})}(\mathbf{x})/\partial {x}_{j}+{f}_{ i}(\mathbf{x})=0$, from equation (12) we see that the spatially varying eigenstrain (4) can be interpreted as a source of an effective body force [8]

$$\begin{equation}{f}_{ i}(\mathbf{x})=-{C}_{ipmn}\frac{\partial {\omega }_{mn}(\mathbf{x})}{\partial {x}_{p}}=-\frac{\partial }{\partial {x}_{p}}{{\Pi}}_{ip}(\mathbf{x}).\end{equation} \tag{ 13 }$$

This interpretation takes the ideal crystalline state of the material as a reference. In this picture, radiation swelling originates from the internal body forces generated by the distribution of defects defined by equation (13). These forces act on the ideal atomic structure of the material, causing it to deform.

A convenient feature of pure elastic strain (9) is that the integral of the trace of the total strain over the volume V of the region occupied by the material, subject to traction-free boundary conditions, equals the sum of relaxation volumes of all the defects in it [8, 43, 50]

$$\begin{equation}{\int }_{V}{{\epsilon}}_{ii}^{(\text{tot})}(\mathbf{x}){\mathrm{d}}^{3}x={\int }_{V}{\omega }_{ii}(\mathbf{x}){\mathrm{d}}^{3}x.\end{equation} \tag{ 14 }$$

The trace of pure elastic strain vanishes after integration over the same region

$$\begin{equation}{\int }_{V}{{\epsilon}}_{ii}(\mathbf{x}){\mathrm{d}}^{3}x=0.\end{equation} \tag{ 15 }$$

This condition shows that the trace of pure elastic strain _ii (x) either vanishes identically everywhere inside the material, or changes sign, varying from compression to tension. A detailed proof, given in appendix A, shows that _ij (x) in fact satisfies an even stronger condition

$$\begin{equation}{\int }_{V}{{\epsilon}}_{ij}(\mathbf{x}){\mathrm{d}}^{3}x=0.\end{equation} \tag{ 16 }$$

The latter is particularly significant since it shows that it is not only the elastic strain that averages to zero over the volume of the material. The elastic stress, which is a linear combination of elements of the elastic strain tensor, must also vanish after the integration over the volume of the material, provided that the traction-free boundary conditions are satisfied; namely

$$\begin{equation}{\int }_{V}{\sigma }_{ij}(\mathbf{x}){\mathrm{d}}^{3}x=0.\end{equation} \tag{ 17 }$$

Below, we show that this condition of vanishing mean stress proves highly significant in applications.

It is the total local strain ${{\epsilon}}_{ij}^{(\text{tot})}(\mathbf{x})$ that is observed in x-ray diffraction experiments. For example, a negative value of total strain, signifying lattice contraction, is observed in an elastically unstrained material containing high spatially homogeneous concentration of vacancies [51]. On the other hand, large positive strain associated with the accumulation of self-interstitial atom defects is observed in thin near-surface layers of tungsten implanted with self-ions to a relatively low dose [25]. The integral of the trace of the total strain ${{\epsilon}}_{ii}^{(\text{tot})}(\mathbf{x})$ over the geometric volume of the material determines the magnitude of dimensional changes and the amount of volumetric swelling occurring following exposure to radiation, see appendix B for more detail.

We highlight the macroscopic nature of distinction between the total and pure elastic components of strain. At a microscopic level, where we are interested in the local deformation of the lattice in the vicinity of a particular defect, there is no need to separate the total and elastic components of strain. However, at the macroscale, where the density of relaxation volumes of defects (4) becomes a meaningful macroscopic quantity, separating the total and elastic components of strain is fully warranted and in fact necessary.

In what follows, we explore several examples illustrating how elastic strain and stress, (9) and (10), vary in reactor components exposed to irradiation. We find elastic strain and stress fields in a rectangular-shape component where the distribution of defects is either spatially homogeneous or spatially heterogeneous. This geometry is also commonly used in ion implantation experiments.

We also investigate the analytically tractable cases of components with cylindrical and spherical geometries exposed to symmetric sources of radiation. Pressurised cylindrical pipes are used in neutron irradiation tests probing how materials respond to exposure to neutrons in the presence of external stress, generated by the internal gas pressure in a pipe. A spherical shell provides a convenient generic case for exploring radiation effects in materials since the displacement, strain, and stress fields turn out to depend only on one, radial, variable. This enables finding exact analytical solutions suitable for comparison with numerical FEM analysis. Abaqus 2021 computer programme was used for this numerical work. The case studies given below enable exploring the effect of spatial variation of defect densities, stemming from the variation of neutron flux, on the magnitude of elastic strain and stress developing in reactor components following their exposure to irradiation.

We find that in the absence of constraints imposed by other elements of the reactor structure, a spatially homogeneous distribution of defects gives rise to volumetric swelling but produces no elastic stress in a component. On the other hand, a strongly spatially heterogeneous distribution of defects, stemming from a strongly spatially varying radiation exposure, may give rise to stress concentrations of significant magnitude, of both compressive and tensile character.

In this section, we provide the general formulae describing the elastic field of point defects produced in a material by exposure to irradiation. We assume the validity of isotropic elasticity approximation, and also that the defects adopt statistically random orientations [8], resulting is that the individual relaxation volume tensors acquire the form

$$\begin{equation}{{\Omega}}_{mn}^{(a)}=\frac{1}{3}{{\Omega}}^{(a)}{\delta }_{mn},\end{equation} \tag{ 18 }$$

where δ_mn is the Dirac delta-symbol, satisfying the condition δ_nn = 3. As a result, the spatially-dependent eigenstrain averaged over all orientations is also diagonal with respect to indices m and n, namely

$$\begin{equation}{{\epsilon}}_{mn}^{\ast }(\mathbf{x})={\omega }_{mn}(\mathbf{x})=\frac{1}{3}\omega (\mathbf{x}){\delta }_{mn},\end{equation} \tag{ 19 }$$

where ω(x) is the density of relaxation volumes of defects [8]. Multiplying (19) by the tensor of elastic constants of the material, which in the isotropic approximation has the form [39, 45]

$$\begin{equation}{C}_{ijkl}=\mu \frac{2\nu }{1-2\nu }{\delta }_{ij}{\delta }_{kl}+\mu \left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right),\end{equation} \tag{ 20 }$$

where ν is the Poisson ratio and μ is the shear modulus, we find the eigenstress [39] generated by the defects

$$\begin{equation}{\sigma }_{ij}^{\ast }(\mathbf{x})={C}_{ijkl}{{\epsilon}}_{kl}^{\ast }(\mathbf{x})=B\omega (\mathbf{x}){\delta }_{ij}.\end{equation} \tag{ 21 }$$

In the above equation, B is the bulk modulus of the material [8, 52]

$$\begin{equation}B=\frac{2\mu }{3}\left(\frac{1+\nu }{1-2\nu }\right).\end{equation} \tag{ 22 }$$

The field of displacements generated by an isotropic point defect with relaxation volume Ω situated at point x' is [50, 53]

$$\begin{align}\hfill {u}_{i}(\mathbf{x})& =\frac{{\Omega}}{12\pi }\left(\frac{1+\nu }{1-\nu }\right)\frac{{x}_{i}-{x}_{i}^{\prime }}{\vert \mathbf{x}-{\mathbf{x}}^{\prime }{\vert }^{3}}\hfill \\ \hfill & =-\frac{{\Omega}}{12\pi }\left(\frac{1+\nu }{1-\nu }\right)\frac{\partial }{\partial {x}_{i}}\frac{1}{\vert \mathbf{x}-{\mathbf{x}}^{\prime }\vert },\hfill \end{align} \tag{ 23 }$$

which follows from equation (6) taken in the isotropic elasticity approximation, with the density of relaxation volume tensors given by ${\omega }_{mn}(\mathbf{x})=\frac{1}{3}{\Omega}{\delta }_{mn}\delta (\mathbf{x}-{\mathbf{x}}^{\prime })$.

Kossevich [50] gives a detailed analysis of the above formula for the field of displacements produced by an isotropic point defect and discusses the range of its validity. In particular, [50] highlights the fact that the field of displacements given by equation (23) does not take into account the boundary conditions at surfaces, and hence an attempt to evaluate the relaxation volume of a defect directly from (23), neglecting the boundary conditions, produces an unexpectedly wrong result.

Similarly, the field of displacements generated directly by a distribution of point defects and not including the effect of boundary conditions, is

$$\begin{equation}{u}_{i}(\mathbf{x})=\frac{1}{12\pi }\left(\frac{1+\nu }{1-\nu }\right){\int }_{V}\omega ({\mathbf{x}}^{\prime })\frac{{x}_{i}-{x}_{i}^{\prime }}{\vert \mathbf{x}-{\mathbf{x}}^{\prime }{\vert }^{3}}\enspace {\mathrm{d}}^{3}{x}^{\prime }.\end{equation} \tag{ 24 }$$

The total strain can now be computed using equation (11). Comparing equations (24) and (11), we see that the trace of the total strain tensor ${{\epsilon}}_{ii}^{(\text{tot})}(\mathbf{x})$, computed directly from (24), equals

$$\begin{equation}{{\epsilon}}_{ii}^{(\text{tot})}(\mathbf{x})=\frac{1}{3}\left(\frac{1+\nu }{1-\nu }\right)\omega (\mathbf{x}).\end{equation} \tag{ 25 }$$

The fact that the right-hand side of the above equation is not equal to ω(x) illustrates the fact that the total change of volume of the material due to the presence of defects is only partially associated with the direct sources of strain, i.e. the local deformations of the lattice produced by the defects, and that partially it arises from the elastic fields formed as a result of application of boundary conditions [43, 50].

A common example of eigenstrain is thermal strain. Equation (1.1) by Mura [39] states that in the isotropic approximation ${{\epsilon}}_{ij}^{\ast }(\mathbf{x})={\delta }_{ij}\alpha T(\mathbf{x})$, where α is the linear thermal expansion coefficient and T(x) is the temperature, provided that the body is in its reference state at T = 0. Applying the condition of mechanical equilibrium, we see that the above thermal strain produces an effective body force [49]

$$\begin{equation}{f}_{ i}=-3B\alpha \frac{\partial T}{\partial {x}_{i}}.\end{equation} \tag{ 26 }$$

Inserting (18) and (20) into (13) and comparing it with (26) shows a way of evaluating stress and strain fields through the use of a fictitious thermal expansion term in FEM, defined by α = 1/3 and T(x) = ω(x). A similar approach was followed by Nguyen et al [54] and Leide et al [55].

Concluding this section, we note that the isotropic formula for the volume tensor of a defect (18) used in a phenomenological treatment of swelling [8, 54–56] does not represent a generally valid approximation. There are documented cases showing that externally applied stress, or the stress developing in a material under irradiation, can polarise defects and generate eigenstrain not diagonal with respect to its Cartesian components [25, 40, 41, 57]. In this study, we do not explore the effects of elastic polarization of defects, but note that equations (9)–(13) enable the evaluation of elastic deformation, strain and stress for an arbitrary anisotropic distribution of eigenstrain produced by the defects generated in materials by irradiation.

Ion irradiation is used as a cost- and time-effective surrogate for neutron irradiation when studying changes in materials properties in a reactor environment. High energy ions have a shorter penetration depth compared to neutrons, and the damaged layer in an ion-irradiated sample is typically only a few microns in thickness. The implantation profile depends on the energy of the ions and can be estimated using the SRIM [58, 59] Monte Carlo simulation code. Figure 1 shows the displacement damage profile (in displacements per atom, or dpa) for W ions implanted into W. The data points shown were obtained from SRIM via the so-called Kinchin–Pease 'quick' calculation (as recommended in [60, 61]). For ion energies of 10, 20 and 50 MeV, the implantation of 10 000 W ions were simulated and the resulting vacancies/ion as a function of depth into the sample were used to calculate the dpa for an incident ion fluence of 10¹⁴ ions cm⁻², which is typical of low-dpa ion-irradiation experiments (see, for example, table S1 in the supplemental material of [25]). Note that, for these SRIM simulations, the ASTM standard threshold displacement energy of 90 eV (see, e.g. [62], table II) was assumed for W.

Zoom In Zoom Out Reset image size

Because of the short ion implantation range, samples for ion-irradiation studies are commonly chosen as thin foils, with micro-mechanical tests subsequently applied to characterize the mechanical properties of the highly stressed damaged layer. For this purpose, it is desirable to be able to fully quantify the relation between the residual stress, or eigenstress, and the elastic fields in the sample. In the limit where the foil thickness is small relative to its width, we are able to find analytical solutions for the elastic stress and strain fields resulting from the accumulation of defects in the foil.

The general solution for the field of displacements is obtained by solving for the condition of mechanical equilibrium (12), which is a linear inhomogeneous partial differential equation. In what follows, we solve the problem using a standard approach by first finding a particular solution of equation (24), and then by finding a solution to the homogeneous problem, which follows from the condition of mechanical equilibrium (12) for ω_mn (x) = 0. The full solution is given by the sum of the particular and homogeneous solutions, where the unknown constants entering the homogeneous solution are identified from the boundary conditions.

Consider a foil with its planar surface parallel to the xy-plane, irradiated by ions coming from the z-direction, producing a relaxation volume density profile ω(z). As the foil width in x- and y-directions is considered to be large relative to its thickness, and the foil is attached to a substrate such that it is unable to buckle, we assume that the strain and stress fields are independent of x and y coordinates. Taking the ensemble-averaged density of relaxation volume tensors (5) as diagonal with respect to the Cartesian indices [8], we write

$$\begin{equation}{\omega }_{mn}(\mathbf{x})=\frac{1}{3}\omega (z){\delta }_{mn}.\end{equation} \tag{ 27 }$$

The field of displacements is defined on the interval 0 ⩽ z ⩽ h, where h is the thickness of the foil. The particular solution, excluding the effect of boundary conditions, follows immediately from (25) as

$$\begin{equation}\frac{\partial }{\partial z}{u}_{z}(z)=\frac{1}{3}\left(\frac{1+\nu }{1-\nu }\right)\omega (z),\end{equation} \tag{ 28 }$$

where we noted that the in-plane displacements vanish in the thin film limit, as can be confirmed by explicitly integrating (24). The field of displacements can now be found by integrating (28) with respect to z, namely

$$\begin{equation}{u}_{z}(z)=\frac{1}{3}\left(\frac{1+\nu }{1-\nu }\right){\int }_{0}^{z}\omega (\zeta )\mathrm{d}\zeta .\end{equation} \tag{ 29 }$$

We now need to find a solution to the homogeneous problem. Since the strain field does not depend on x or y, the field of displacements has the form

$$\begin{equation}\mathbf{u}(\mathbf{x})=cx{\mathbf{e}}_{x}+cy{\mathbf{e}}_{y}+{u}_{z}(z){\mathbf{e}}_{z},\end{equation} \tag{ 30 }$$

where e_x , e_y and e_z are the Cartesian unit vectors, and c is a constant. Substituting (30) into the homogeneous condition of mechanical equilibrium, we arrive at an ordinary differential equation

$$\begin{equation}\frac{{\mathrm{d}}^{2}{u}_{z}(z)}{\mathrm{d}{z}^{2}}=0.\end{equation} \tag{ 31 }$$

This equation has general solutions of the form

$$\begin{equation}{u}_{z}(z)=a+bz,\end{equation} \tag{ 32 }$$

where we are free to set a ≡ 0 as uniform displacements do not affect elastic fields.

A general solution for u(x), valid in the interval 0 ⩽ z ⩽ h, is given by the sum of the homogeneous solution (32) and the particular solution (29), namely

$$\begin{equation}\begin{aligned}\hfill {u}_{x}(x)& =cx\hfill \\ \hfill {u}_{y}(y)& =cy\hfill \\ \hfill {u}_{z}(z)& =bz+\frac{1}{3}\left(\frac{1+\nu }{1-\nu }\right){\int }_{0}^{z}\omega (\zeta )\mathrm{d}\zeta .\hfill \end{aligned}\end{equation} \tag{ 33 }$$

Combining equations (9), (10), (21) with the above expression (33), we can express the non-vanishing components of the elastic stress tensor as

$$\begin{equation*}\begin{aligned}\hfill {\sigma }_{xx}(z)={\sigma }_{yy}(z)& =\frac{2\mu (c+b\nu )}{1-2\nu }-\frac{2\mu }{3}\left(\frac{1+\nu }{1-\nu }\right)\omega (z),\hfill \\ \hfill {\sigma }_{zz}(z)& =\frac{2\mu b(1-\nu )}{1-2\nu }+\frac{4\mu c\nu }{1-2\nu }.\hfill \end{aligned}\end{equation*}$$

All the other components of tensor σ_ij (z) are equal to zero.

Constants b and c depend on the conditions at the top and bottom surfaces of the foil. In the context of an ion-irradiation experiment, the surface of the foil at z = h, facing the ion flux, is free of tractions, therefore the first condition is σ_zz (h) = 0. For the second condition, we consider first the case where the foil is able to relax freely in the xy-plane. While it is not possible to derive a solution where the foil sides are completely free of tractions, as that would require buckling of the foil, we can weaken the condition and instead assert that the average traction at the sides of the foil vanishes, which corresponds to a sliding condition. For example, the average traction over the foil boundaries at x → ±∞ is

$$\begin{equation}\frac{1}{h}{\int }_{0}^{h}{t}_{i}\enspace \mathrm{d}z=\frac{1}{h}{\int }_{0}^{h}{\sigma }_{ij}{n}_{j}\enspace \mathrm{d}z=\frac{1}{h}{\int }_{0}^{h}{\sigma }_{xx}(z)\mathrm{d}z=0,\end{equation} \tag{ 34 }$$

where n = e_x is the outward surface normal vector to the side of the foil. Solving for b and c such that conditions σ_zz (h) = 0 and (34) are met, we arrive at explicit analytical expressions for the field of displacements

$$\begin{equation}\begin{aligned}\hfill {u}_{x}(x)& =\frac{1}{3}\bar{\omega }x\hfill \\ \hfill {u}_{y}(y)& =\frac{1}{3}\bar{\omega }y\hfill \\ \hfill {u}_{z}(z)& =-\frac{2\nu }{3(1-\nu )}\bar{\omega }z+\frac{1}{3}\left(\frac{1+\nu }{1-\nu }\right){\int }_{0}^{z}\omega (\zeta )\mathrm{d}\zeta ,\hfill \end{aligned}\end{equation} \tag{ 35 }$$

where $\bar{\omega }$ is the mean value of function ω(z) over the thickness of the foil

$$\begin{equation}\bar{\omega }=\frac{1}{h}{\int }_{0}^{h}\omega (\zeta )\mathrm{d}\zeta .\end{equation} \tag{ 36 }$$

The only elastic stress components that do not vanish are the in-plane stresses, which equal

$$\begin{equation}{\sigma }_{xx}(z)={\sigma }_{yy}(z)=-\frac{2\mu }{3}\left(\frac{1+\nu }{1-\nu }\right)\left[\omega (z)-\bar{\omega }\right].\end{equation} \tag{ 37 }$$

Let us now consider the case where the sides of the foil are fixed such that no in-plane displacements occur. The displacement field of the constrained foil ${u}_{i}^{C}(\mathbf{x})$ with the condition ${u}_{x}^{C}={u}_{y}^{C}=0$ is obtained from the sliding foil solution (35) by setting $\bar{\omega }=0$, irrespective of its actual value, as then the in-plane displacements vanish. The field of vertical displacements is now given by

$$\begin{equation}{u}_{z}^{C}(z)=\frac{1}{3}\left(\frac{1+\nu }{1-\nu }\right){\int }_{0}^{z}\omega (\zeta )\mathrm{d}\zeta .\end{equation} \tag{ 38 }$$

The in-plane components of elastic stress are

$$\begin{equation}{\sigma }_{xx}^{C}(z)={\sigma }_{yy}^{C}(z)=-\frac{2\mu }{3}\left(\frac{1+\nu }{1-\nu }\right)\omega (z).\end{equation} \tag{ 39 }$$

Incidentally, the same solution is found by taking the limit h → ∞ in the solution for a sliding foil (35), provided that the damage profile is localised such that ${\mathrm{lim}}_{h\to \infty }\enspace \bar{\omega }=0$. This demonstrates that the in-plane stresses in the solution for a sliding foil approach those of the constrained foil as the foil thickness increases.

Solutions found for the constrained and unconstrained cases show that swelling induces stress only in the plane of the foil, and there is no stress in the direction normal to the foil surface, either at the surface or in the bulk of the foil. The in-plane stress can be fairly high, for example swelling of 1% in a constrained tungsten foil produces an in-plane stress approaching −1.9 GPa.

To verify the analytical solutions given above, we compare them to numerical FEM simulations. We considered a tungsten foil (μ = 160 GPa and ν = 0.28) with dimensions of 80 μm in the xy-plane and height of h = 10 μm along z. Two types of boundary conditions were considered at the bottom surface of the foil interfacing the substrate: the foil is either completely constrained (u_x (0) = u_y (0) = u_z (0) = 0) or its displacements are constrained only in the z direction (u_z (0) = 0). The FEM solutions involved approximately 1.20 × 10⁶ hexahedral linear elements. The free surface of the foil before irradiation is at z = h, and the assumed profile ω(z) is shown in figure 1 as a function of the distance from the top surface, or h − z. Figure 2, where the density of relaxation volumes corresponding to 20 MeV ion irradiation was used as input, shows that the two analytical solutions derived in this section approach the two types of boundary conditions considered for the FEM numerical solutions.

Zoom In Zoom Out Reset image size

It is desirable to quantify the magnitude of elastic stress that develops in a foil following ion irradiation. Indeed, the ion irradiation data can help extrapolate experimental observations to the neutron irradiated case, where the distribution of relaxation volume density ω(x) varies on the scale of tens of centimetres [23]. Figure 3 shows how the in-plane components of stress vary depending on the energy of the incident ions, which produce damage at greater depth with higher energy. The boundary condition at the bottom surface is taken to be the sliding one, allowing for expansion of the film in the xy-plane. The σ_zz component of stress was found to be vanishingly small in all cases (lower than 1 MPa), confirming the results of analytical treatment. The in-plane compressive stress components are significant, and exceed −1 GPa in the irradiated layer even if the foil is able to relax and expand sideways.

Zoom In Zoom Out Reset image size

The development of compressive or dilatational stresses in the GPa range has been observed in a variety of materials exposed to ion irradiation, including ceramic oxides (−2.7 GPa in Gd₂TiZrO₇ at 0.2 dpa) [63], SiC and SiO₂ [64], Cr [65], Au [66] and W [25] thin films. Davis [67] estimated the compressive stress in thin films growing on a substrate, assuming that a steady state is attained when the rate of implantation of new ions is balanced by the rate with which stress is relaxed by atoms escaping to the surface. Wolfer and Garner proposed instead that the compressive stress depends on the competition between swelling and creep [68]. These examples show that taking into account the effects of stress relaxation—whether by thermal or irradiation creep—is required before comparing theoretical predictions with experimental results.

A detailed study of strain and stress in ion-irradiated tungsten [25, 26] shows that the compressive stress in the ion-irradiated layer initially reaches ∼320 MPa [69] at 0.3 dpa before changing sign and equilibrating, at an approximately similar level in magnitude but dilatational in character, in the high dose limit. This finding agrees with the earlier observations by Misra et al [65] and illustrates the highly non-linear nature of the stress relaxation phenomenon. It appears that the radiation-driven relaxation still does not fully eliminate the stress produced by the defects themselves, resulting in the residual stress in the several 100 MPa range [16, 25].

In this section we investigate the distribution of elastic strain and stress in a sample of cylindrical geometry. Cylindrical tubes are often used in the context of neutron irradiation experiments not only because of the natural similarity with the geometry of fission fuel cladding but also because a sample of cylindrical geometry offers a convenient way of testing the combined effect of externally applied stress and irradiation. In practice, the external stress is produced by the pressure of gas inside the cylinder, whereas the exposure to irradiation is often external.

Consider a cylindrical tube with the inner radius ρ₁ and outer radius ρ₂, and define a cylindrical coordinate system (ρ, φ, z) coaxial with the tube. Pressure of gas P applied to the inner surface is assumed to be constant, and the density of relaxation volumes tensor is a function of only the radial coordinate ρ, namely

$$\begin{equation}{\omega }_{mn}(\mathbf{x})=\frac{\omega (\rho )}{3}{\delta }_{mn}.\end{equation} \tag{ 40 }$$

Since neither P nor ω(ρ) depend on φ or z, the field of displacements has the form

$$\begin{equation}\mathbf{u}(\mathbf{x})={u}_{\rho }(\rho ){\boldsymbol{e}}_{\boldsymbol{\rho }}+cz{\boldsymbol{e}}_{z},\end{equation} \tag{ 41 }$$

where c is a constant and e _ρ and e_z are the unit vectors in the radial and z directions, respectively. The choice of the above simple form of the field of displacements stems from the fact that u_φ = 0 from symmetry, and the axial strain must be translationally invariant with respect to z, implying that ∂u_z /∂z = c. For a pipe of finite length, this translational invariance only applies to the central portion of the pipe, far from its two extremities. Below, we explore two different boundary conditions along z, namely that the pipe is either fully unconstrained or fully constrained. We find that the irradiation-induced swelling produces different patterns of stress in these two limits.

If the tube is constrained in the z direction, then c = 0. If the distribution of ω_mn (x) is entirely spatially uniform, the components of stress can be found by solving the homogeneous equation ∂σ_ij /∂x_j = 0. Evaluating the divergence of stress in the cylindrical system of coordinates in the isotropic elasticity approximation, the field of radial displacements can be found by solving the ordinary differential equation

$$\begin{equation}\rho \frac{{\mathrm{d}}^{2}{u}_{\rho }}{\mathrm{d}{\rho }^{2}}+\frac{\mathrm{d}{u}_{\rho }}{\mathrm{d}\rho }-\frac{{u}_{\rho }}{\rho }=0.\end{equation} \tag{ 42 }$$

A general solution of this homogeneous equation has the form [49]

$$\begin{equation}{u}_{\rho }(\rho )=a\rho +\frac{b}{\rho }.\end{equation} \tag{ 43 }$$

To find a particular solution describing the effect of relaxation volumes of defects given by equation (40), we note that the divergence of the displacement field generated by ω( ρ ) is given by equation (25), namely

$$\begin{equation}\frac{\partial }{\partial {x}_{i}}{u}_{i}(\boldsymbol{\rho })={{\epsilon}}_{ii}^{(\text{tot})}(\mathbf{x})=\frac{1}{3}\left(\frac{1+\nu }{1-\nu }\right)\omega (\boldsymbol{\rho }).\end{equation} \tag{ 44 }$$

Using the divergence theorem, and noting that the field of displacements depends only on the radial variable, equation (44) gives the particular solution in the form

$$\begin{equation}{u}_{\rho }(\rho )=\frac{1}{3\rho }\left(\frac{1+\nu }{1-\nu }\right)\underset{{\rho }_{1}}{\overset{\rho }{\int }}R\omega (R)\mathrm{d}R.\end{equation} \tag{ 45 }$$

The general solution, valid for ρ₁ < ρ < ρ₂, is given by the sum of the homogeneous solution (43) and the particular solution (45), namely

$$\begin{equation}{u}_{\rho }(\rho )=a\rho +\frac{b}{\rho }+\frac{1}{3\rho }\left(\frac{1+\nu }{1-\nu }\right)\underset{{\rho }_{1}}{\overset{\rho }{\int }}R\omega (R)\mathrm{d}R.\end{equation} \tag{ 46 }$$

Combining equations (9), (10), and (21), the components of the stress tensor can now be expressed as

$$\begin{equation*}\begin{aligned}\hfill {\sigma }_{\rho \rho }(\rho )& =2\mu {{\epsilon}}_{\rho \rho }(\rho )+\frac{2\mu \nu }{1-2\nu }{{\epsilon}}_{ii}(\rho )-B\omega (\rho )\hfill \\ \hfill {\sigma }_{\varphi \varphi }(\rho )& =2\mu {{\epsilon}}_{\varphi \varphi }(\rho )+\frac{2\mu \nu }{1-2\nu }{{\epsilon}}_{ii}(\rho )-B\omega (\rho )\hfill \\ \hfill {\sigma }_{zz}(\rho )& =2\mu {{\epsilon}}_{zz}+\frac{2\mu \nu }{1-2\nu }{{\epsilon}}_{ii}(\rho )-B\omega (\rho ),\hfill \end{aligned}\end{equation*}$$

where _ii = _ρρ + _φφ + _zz is the trace of the elastic strain tensor and _zz is a constant independent of ρ.

In the unconstrained case, we still have to determine the three constants a, b, and c in expression (41), derived from equation (46). Two of the three conditions required for determining the three constants stem from the boundary conditions σ_ρρ (ρ₁) = −P, where P is the pressure applied to the internal surface, and σ_ρρ (ρ₂) = 0. The third condition can be obtained from the condition of mechanical equilibrium. If a body has an external surface S loaded by external traction forces t_i and is arbitrarily sectioned, then the cut defines an internal section Σ with normal vector n_j that divides the external surface into S₁ and S₂. The equilibrium condition applied to the part of the volume enclosed by S₁ and Σ implies that

$$\begin{equation*}{\int }_{{S}_{1}}{t}_{i}\enspace \mathrm{d}{S}_{1}={\int }_{{\Sigma}}{\sigma }_{ij}{n}_{j}\enspace \mathrm{d}{\Sigma}.\end{equation*}$$

Noting that t_i = 0 since the pipe is free of external loads, and taking a section perpendicular to z, we see that

$$\begin{equation*}2\pi {\int }_{{\rho }_{1}}^{{\rho }_{2}}{\sigma }_{zz}(R)\mathrm{d}R=0,\end{equation*}$$

that is, σ_zz must be self-equilibrated. Combining the three boundary conditions, we find the radial and axial displacements

$$\begin{align}\hfill {u}_{\rho }(\rho )& =\frac{P}{2\mu \rho }\frac{\pi {\rho }_{1}^{2}}{A}\left[{\rho }_{2}^{2}+\left(\frac{1-\nu }{1+\nu }\right){\rho }^{2}\right]+\frac{1}{6\rho }\left(\frac{1+\nu }{1-\nu }\right)\hfill \\ \hfill & \quad \times \left[{\rho }_{1}^{2}\bar{\omega }+\left(\frac{1-3\nu }{1+\nu }\right){\rho }^{2}\bar{\omega }+2{\int }_{{\rho }_{1}}^{\rho }R\omega (R)\mathrm{d}R\right]\hfill \\ \hfill {u}_{z}(z)& =\left[\frac{\bar{\omega }}{3}-\frac{P}{\mu }\frac{\pi {\rho }_{1}^{2}}{A}\left(\frac{\nu }{1+\nu }\right)\right]z,\hfill \end{align} \tag{ 47 }$$

where $\bar{\omega }$ is the mean value of function ω(r) over the cross-section area of the cylinder $A=\pi \left({\rho }_{2}^{2}-{\rho }_{1}^{2}\right)$,

$$\begin{equation}\bar{\omega }=\frac{1}{A}{\int }_{S}\omega (\mathbf{x})\mathrm{d}S=\frac{2\pi }{A}{\int }_{{\rho }_{1}}^{{\rho }_{2}}R\omega (R)\mathrm{d}R.\end{equation} \tag{ 48 }$$

Noting that _ρρ = ∂u_ρ /∂ρ, _φφ = u_ρ /ρ, and _zz = ∂u_z /∂z, we find the three non-vanishing components of stress

$$\begin{align}\hfill {\sigma }_{\rho \rho }(\rho )& =\frac{\pi {\rho }_{1}^{2}P}{A}\left(1-\frac{{\rho }_{2}^{2}}{{\rho }^{2}}\right)+\frac{\mu }{3}\left(\frac{1+\nu }{1-\nu }\right)\hfill \\ \hfill & \quad \times \left[\left(1-\frac{{\rho }_{1}^{2}}{{\rho }^{2}}\right)\bar{\omega }-\frac{2}{{\rho }^{2}}{\int }_{{\rho }_{1}}^{\rho }R\omega (R)\mathrm{d}R\right],\hfill \end{align} \tag{ 49 }$$

$$\begin{align}\hfill {\sigma }_{\varphi \varphi }(\rho )& =\frac{\pi {\rho }_{1}^{2}P}{A}\left(1+\frac{{\rho }_{2}^{2}}{{\rho }^{2}}\right)+\frac{\mu }{3}\left(\frac{1+\nu }{1-\nu }\right)\hfill \\ \hfill & \quad \times \left[\left(1+\frac{{\rho }_{1}^{2}}{{\rho }^{2}}\right)\bar{\omega }+\frac{2}{{\rho }^{2}}{\int }_{{\rho }_{1}}^{\rho }R\omega (R)\mathrm{d}R-2\omega (\rho )\right],\hfill \end{align} \tag{ 50 }$$

$$\begin{equation}{\sigma }_{zz}(\rho )=\frac{2\mu }{3}\left(\frac{1+\nu }{1-\nu }\right)\left[\bar{\omega }-\omega (\rho )\right].\end{equation} \tag{ 51 }$$

The above equations show that σ_ρρ (ρ) in (49) vanishes at ρ = ρ₁ and ρ = ρ₂ for an arbitrary distribution ω(ρ). Because of the linear superposition principle, effects of internal pressure and irradiation are completely decoupled. In the absence of irradiation, the expressions for the components of stress readily reduce to the Lamé equations for a thick-walled pressurised cylinder [70]. If swelling, although present, is spatially homogeneous, its contribution to all the stress components vanishes identically, as expected for an unconstrained body subject to traction-free boundary conditions.

If the tube is constrained along z, so that u_z = 0, _zz = 0, then only two conditions need to be specified to determine the constants in equation (46). It is sufficient to prescribe the values of the radial stress at internal and external surfaces, namely σ_ρρ (ρ₁) = −P and σ_ρρ (ρ₂) = 0. In fact, the presence of external forces, required to keep the length of the tube constant, means that σ_zz need not to self-equilibrate. The radial displacement function in the constrained (C) case is

$$\begin{align}\hfill {u}_{\rho }^{C}(\rho )& =\frac{P}{2\mu \rho }\frac{\pi {\rho }_{1}^{2}}{A}\left[{\rho }_{2}^{2}+(1-2\nu ){\rho }^{2}\right]+\frac{1}{6\rho }\left(\frac{1+\nu }{1-\nu }\right)\hfill \\ \hfill & \quad \times \left[{\rho }_{1}^{2}\bar{\omega }+(1-2\nu ){\rho }^{2}\bar{\omega }+2{\int }_{{R}_{1}}^{\rho }R\omega (R)\mathrm{d}R\right].\hfill \end{align} \tag{ 52 }$$

We find that the stresses σ_ρρ (ρ) and σ_φφ (ρ) are not affected by the different boundary conditions and remain the same as in equations (49) and (50). σ_zz (ρ) in the constrained case, on the other hand, has the form

$$\begin{equation}{\sigma }_{zz}^{C}(\rho )=\frac{2\pi \nu {\rho }_{1}^{2}P}{A}+\frac{2\mu }{3}\left(\frac{1+\nu }{1-\nu }\right)\left[\nu \bar{\omega }-\omega (\rho )\right].\end{equation} \tag{ 53 }$$

A comparison of equations (51) and (53) highlights an important point. Spatially homogeneous irradiation generates no σ_zz component of elastic stress if the pipe is unconstrained, but if the extremities are fixed and P = 0, then we have

$$\begin{equation}{\sigma }_{zz}^{C}=-\frac{2\mu }{3}(1+\nu )\bar{\omega }.\end{equation} \tag{ 54 }$$

The above equations shows that the magnitude of elastic stress generated by irradiation is of the order of ω times the shear modulus. As the former can be as large as 1% or more, see [19, 20], internal elastic stress reaching hundreds of MPa can develop in the absence of any internal gas pressure in a pipe with constrained ends, even in the limit where irradiation is spatially homogeneous.

Dimensional stability is an important consideration in the context of design and operation of a nuclear reactor involving many interdependent components. As a result of applied pressure and irradiation, the volume of a pipe can change. This volume change is given by the volume integral of the trace of the total strain. If the pipe is free then, far from its extremities, the volume change of the pipe per unit length can be computed as

$$\begin{align}\hfill \frac{{\Delta}V}{L}& =\frac{1}{L}{\int }_{V}\frac{\partial {u}_{i}}{\partial {x}_{i}}\enspace \mathrm{d}V=2\pi {\int }_{{\rho }_{1}}^{{\rho }_{2}}\left(\frac{\partial (\rho {u}_{\rho })}{\partial \rho }+\rho \frac{\partial {u}_{z}}{\partial z}\right)\mathrm{d}\rho \hfill \\ \hfill & =\frac{P}{\mu }\left(\frac{1-2\nu }{1+\nu }\right)\pi {\rho }_{1}^{2}+A\bar{\omega }.\hfill \end{align} \tag{ 55 }$$

If the pipe is constrained, the change of volume is given by

$$\begin{align}\hfill \frac{{\Delta}{V}^{C}}{L}& =\frac{1}{L}{\int }_{V}\frac{\partial {u}_{i}^{C}}{\partial {x}_{i}}\enspace \mathrm{d}V\hfill \\ \hfill & =\frac{P}{\mu }(1-2\nu )\pi {\rho }_{1}^{2}+\frac{2}{3}(1+\nu )A\bar{\omega }.\hfill \end{align} \tag{ 56 }$$

Equation (55) shows that in the limit where P = 0, the volume change of an unconstrained irradiated pipe equals the total relaxation volume of all the defects in the material. If the pipe is constrained, the volume change is smaller than the total relaxation volume of defects, and the difference can be attributed to the effect of external constraints. Both equations (55) and (56) confirm the intuitive conclusion that in the limit where the material is incompressible (ν = 1/2), the application of external pressure P does not alter the volume of the pipe.

Equations (55) and (56) can also be derived by evaluating, to first order, the change of volume associated with the transformation of a cylinder from its initial undeformed configuration to the final one with the inner radius of ρ₁ + u_ρ (ρ₁), outer radius of ρ₂ + u_ρ (ρ₂) and the length of L + u_z (L).

An FEM model was constructed for a pipe made of ferritic steel (μ = 80 GPa, ν = 0.29), with ρ₁ = 2.5 mm, ρ₁ = 3.0 mm and pressure of P = 50 MPa acting on its inner surface. These parameters are close to the experimental conditions for pressurised steel specimens used in irradiation tests [15]. Function ω(ρ) was selected in such a way so that the density of relaxation volumes varies from about 1% at the outer surface, directly exposed to irradiation, to about 0.5% at the inner surface, see equation (C2). Numerical FEM solutions are illustrated in figure 4. Simulations were performed assuming the effect of pressure only, irradiation only, and the combined effect of both pressure and irradiation. To achieve agreement between numerical FEM results and exact analytical solutions, about 9.40 × 10⁵ hexahedral linear elements were used.

Zoom In Zoom Out Reset image size

If the exposure of a pipe to an external source of irradiation gives rise to a spatially varying defect density, the maximum tensile hoop stress, σ_φφ , occurs at the inner surface of the pipe, which is also where the hoop stress induced by P is maximum. The axial stress, σ_zz , is negative throughout the thickness of the pipe if it is constrained. If the pipe can expand freely along its length, then the axial stress must change sign to satisfy the condition of self-equilibrium. A change in the density of relaxation volumes of 0.5% produces stresses that are comparable with those generated by high internal pressure, with potential implication for the design of experimental tests.

In this section we evaluate elastic stress and strain developing in a spherical shell exposed to irradiation. The advantage offered by the spherical symmetry of the shell is that in the limit where the source of irradiation is also spherically symmetric, the solutions can be found in a closed analytical form. Some of them, related to the total strain and stress, were investigated earlier in reference [8], however the pure elastic components of stress and strain were not evaluated. Given the significance of pure elastic solutions in the context of assessment of structural integrity of power plant components, we present the relevant analysis below.

If the source of irradiation is spherically symmetric, the distribution of defect relaxation volumes depends only on the distance r to the centre of the shell. It 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 density of relaxation volume tensors now has the form

$$\begin{equation}{\omega }_{mn}(\mathbf{x})=\frac{1}{3}\omega (r){\delta }_{mn},\end{equation} \tag{ 57 }$$

and the field of displacements is defined on the interval R₁ ⩽ r ⩽ R₂, where R₁ is the inner radius of the shell and R₂ is its outer radius. This field is radially-symmetric and the vector of displacements can be written as u(r) = u_r (r)n, where n = r/r. The divergence of the field of atomic displacements has the form (25)

$$\begin{equation}\frac{\partial }{\partial r}u(r)=\frac{1}{3}\left(\frac{1+\nu }{1-\nu }\right)\omega (r).\end{equation} \tag{ 58 }$$

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

$$\begin{equation}4\pi {r}^{2}{u}_{r}(r)=\frac{4\pi }{3}\left(\frac{1+\nu }{1-\nu }\right){\int }_{{R}_{1}}^{r}{R}^{2}\omega (R)\mathrm{d}R,\end{equation} \tag{ 59 }$$

for R₁ ⩽ r ⩽ R₂.

Noting that in the absence of radiation defects the divergence of u(x) is a harmonic function of coordinates, we write the field of displacements as a sum of the particular solution of heterogeneous equation (58) and a general solution of the corresponding homogeneous equation, namely

$$\begin{equation}{u}_{r}(r)=ar+\frac{b}{{r}^{2}}+\frac{1}{3{r}^{2}}\left(\frac{1+\nu }{1-\nu }\right){\int }_{{R}_{1}}^{r}{R}^{2}\omega (R)\mathrm{d}R.\end{equation} \tag{ 60 }$$

Here a and b are constants that need to be determined from the 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 σ_ij (r)n_j = 0 at r = R₁ and r = R₂. The total strain can be found by differentiating (60). It has the form

$$\begin{align}\hfill {{\epsilon}}_{ij}^{(\text{tot})}(r)& =a{\delta }_{ij}+\frac{b}{{r}^{3}}\left({\delta }_{ij}-3{n}_{i}{n}_{j}\right)\hfill \\ \hfill & \quad +\frac{1}{3{r}^{3}}\left(\frac{1+\nu }{1-\nu }\right)\left({\delta }_{ij}-3{n}_{i}{n}_{j}\right){\int }_{{R}_{1}}^{r}{R}^{2}\omega (R)\mathrm{d}R\hfill \\ \hfill & \quad +\frac{1}{3}\left(\frac{1+\nu }{1-\nu }\right){n}_{i}{n}_{j}\omega (r).\hfill \end{align} \tag{ 61 }$$

The find the pure elastic part of strain, we need to subtract the eigenstrain of defects $\frac{1}{3}\omega (r){\delta }_{ij}$ from the above expression, resulting in

$$\begin{align}\hfill {{\epsilon}}_{ij}(r)& =a{\delta }_{ij}+\frac{b}{{r}^{3}}\left({\delta }_{ij}-3{n}_{i}{n}_{j}\right)\hfill \\ \hfill & \quad +\frac{1}{3{r}^{3}}\left(\frac{1+\nu }{1-\nu }\right)\left({\delta }_{ij}-3{n}_{i}{n}_{j}\right){\int }_{{R}_{1}}^{r}{R}^{2}\omega (R)\mathrm{d}R\hfill \\ \hfill & \quad +\frac{1}{3}\left[\left(\frac{1+\nu }{1-\nu }\right){n}_{i}{n}_{j}-{\delta }_{ij}\right]\omega (r).\hfill \end{align} \tag{ 62 }$$

To find the pure elastic stress, we multiply _ij (r) by the forth-rank tensor of elastic constants (20). The resulting expression for elastic stress is

$$\begin{align}\hfill {\sigma }_{ij}(r)& =2\mu a\left(\frac{1+\nu }{1-2\nu }\right){\delta }_{ij}+\frac{2\mu b}{{r}^{3}}({\delta }_{ij}-3{n}_{i}{n}_{j})\hfill \\ \hfill & \quad +\frac{2\mu }{3{r}^{3}}\left(\frac{1+\nu }{1-\nu }\right)\left({\delta }_{ij}-3{n}_{i}{n}_{j}\right){\int }_{{R}_{1}}^{r}{R}^{2}\omega (R)\mathrm{d}R\hfill \\ \hfill & \quad +\frac{\mu }{3}\left(\frac{1+\nu }{1-\nu }\right)\omega (r)\left(\frac{2\nu }{1-2\nu }{\delta }_{ij}+2{n}_{i}{n}_{j}\right)\hfill \\ \hfill & \quad -\frac{2\mu }{3}\left(\frac{1+\nu }{1-2\nu }\right)\omega (r){\delta }_{ij}.\hfill \end{align} \tag{ 63 }$$

Projecting elastic stress onto the radial unit vector n = r/r, we find

$$\begin{align}\hfill {\sigma }_{ij}(r){n}_{j}& =2\mu a\left(\frac{1+\nu }{1-2\nu }\right){n}_{i}-\frac{4\mu b}{{r}^{3}}{n}_{i}\hfill \\ \hfill & \quad -\frac{4\mu }{3{r}^{3}}\left(\frac{1+\nu }{1-\nu }\right){n}_{i}{\int }_{{R}_{1}}^{r}{R}^{2}\omega (R)\mathrm{d}R.\hfill \end{align} \tag{ 64 }$$

Applying the traction-free boundary conditions σ_ij (r)n_j = 0 at r = R₁ and r = R₂ to (63) and (64), we find parameters a and b, namely

$$\begin{equation}a=\frac{2}{9}\left(\frac{1-2\nu }{1-\nu }\right)\bar{\omega },\end{equation} \tag{ 65 }$$

and

$$\begin{equation}b=\frac{1}{9}\left(\frac{1+\nu }{1-\nu }\right)\bar{\omega },\end{equation} \tag{ 66 }$$

where $\bar{\omega }$ is the mean density of relaxation volumes of all the defects accumulated in the shell

$$\begin{equation}\bar{\omega }=\frac{1}{V}{\int }_{V}\omega (\mathbf{x})\mathrm{d}V=\frac{4\pi }{V}{\int }_{{R}_{1}}^{{R}_{2}}{R}^{2}\omega (R)\mathrm{d}R,\end{equation} \tag{ 67 }$$

and V is the geometric volume of the shell

$$\begin{equation}V=\frac{4\pi }{3}\left({R}_{2}^{3}-{R}_{1}^{3}\right).\end{equation} \tag{ 68 }$$

The above formulae remain valid irrespectively of whether or not ω(r) vanishes at surfaces r = R₁ and r = R₂.

The total macroscopic change of volume resulting from the accumulation of defects in the shell equals the integral of the trace of the total strain tensor (14) over the volume of the component

$$\begin{align}\hfill {\Delta}V& ={\int }_{V}{{\epsilon}}_{ii}(\mathbf{r})\mathrm{d}V={\int }_{V}\left[\frac{1}{3}\left(\frac{1+\nu }{1-\nu }\right)\omega (\mathbf{r})+3a\right]\mathrm{d}V\hfill \\ \hfill & ={\int }_{{R}_{1}}^{{R}_{2}}\left[\frac{1}{3}\left(\frac{1+\nu }{1-\nu }\right)\omega (R)+3a\right]4\pi {R}^{2}\enspace \mathrm{d}R,\hfill \end{align} \tag{ 69 }$$

where we noted that δ_ii = 3 and n_i n_i = 1. The integral equals ${{\Omega}}_{\mathrm{t}\mathrm{o}\mathrm{t}}=V\bar{\omega }$, in agreement with equation (14).

Substituting (65) and (66) into (60), we find the radial displacements of the inner and outer surfaces of the shell

$$\begin{equation}\begin{aligned}\hfill {u}_{r}({R}_{1})& =\frac{1}{3}{R}_{1}\bar{\omega }\hfill \\ \hfill {u}_{r}({R}_{2})& =\frac{1}{3}{R}_{2}\bar{\omega }.\hfill \end{aligned}\end{equation} \tag{ 70 }$$

These displacements satisfy the condition

$$\begin{equation}4\pi {R}_{2}^{2}{u}_{r}({R}_{2})-4\pi {R}_{1}^{2}{u}_{r}({R}_{1})=V\bar{\omega }={{\Omega}}_{\mathrm{t}\mathrm{o}\mathrm{t}},\end{equation} \tag{ 71 }$$

which provides an alternative way of evaluating the total volumetric swelling of the shell [8]. Figure 5(a) graphically exemplifies equation (14) and is in full agreement with (71), for two functions ω(r) that have the same $\bar{\omega }$.

Zoom In Zoom Out Reset image size

To assess the performance of a component under operating conditions where radiation defects accumulate in the bulk of its structure, it is convenient to use the spherical system of coordinates and project the elastic stress tensor (63) onto the three orthogonal unit vectors (e_r , e_θ , e_ϕ ), corresponding to the radial and angular degrees of freedom and related to the various modes of deformation of the shell. These unit vectors are

$$\begin{equation}\begin{aligned}\hfill {\mathbf{e}}_{r}& ={\mathbf{e}}_{x}\enspace \mathrm{sin}\enspace \theta \enspace \mathrm{cos}\enspace \phi +{\mathbf{e}}_{y}\enspace \mathrm{sin}\enspace \theta \mathrm{sin}\enspace \phi +{\mathbf{e}}_{z}\enspace \mathrm{cos}\enspace \theta \hfill \\ \hfill {\mathbf{e}}_{\theta }& ={\mathbf{e}}_{x}\enspace \mathrm{cos}\enspace \theta \enspace \mathrm{cos}\enspace \phi +{\mathbf{e}}_{y}\enspace \mathrm{cos}\enspace \theta \enspace \mathrm{sin}\enspace \phi -{\mathbf{e}}_{z}\enspace \mathrm{sin}\enspace \theta \hfill \\ \hfill {\mathbf{e}}_{\phi }& =-{\mathbf{e}}_{x}\enspace \mathrm{sin}\enspace \phi +{\mathbf{e}}_{y}\enspace \mathrm{cos}\enspace \phi ,\hfill \end{aligned}\end{equation} \tag{ 72 }$$

where θ and ϕ are the polar and azimuthal angles of the spherical system of coordinates. Since e_r is the same as n, we have (n ⋅ e_θ ) = (n ⋅ e_ϕ ) = 0.

The radial diagonal element of the elastic stress tensor is

$$\begin{align}\hfill {\sigma }_{rr}(r)& =\frac{4\mu }{9}\left(\frac{1+\nu }{1-\nu }\right)\left[\bar{\omega }\left(1-\frac{{R}_{1}^{3}}{{r}^{3}}\right)\right]\hfill \\ \hfill & \quad -\frac{4\mu }{9}\left(\frac{1+\nu }{1-\nu }\right)\left[\frac{3}{{r}^{3}}{\int }_{{R}_{1}}^{r}{R}^{2}\omega (R)\mathrm{d}R\right].\hfill \end{align} \tag{ 73 }$$

A direct calculation shows that in the limit where function ω(r) does not depend on r, ω(r) = const, the radial component of elastic stress vanishes identically everywhere in the volume of the shell, σ_rr (r) = 0 for all R₁ < r < R₂.

The circumferential (hoop) components σ_θθ (r) and σ_ϕϕ (r) of the elastic stress tensor are

$$\begin{align}\hfill {\sigma }_{\theta \theta }(r)& =\frac{4\mu }{9}\left(\frac{1+\nu }{1-\nu }\right)\left[\bar{\omega }\left(1+\frac{{R}_{1}^{3}}{2{r}^{3}}\right)\right]\hfill \\ \hfill & \quad +\frac{2\mu }{3}\left(\frac{1+\nu }{1-\nu }\right)\left[\frac{1}{{r}^{3}}\underset{{R}_{1}}{\overset{r}{\int }}{R}^{2}\omega (R)\mathrm{d}R\right]\hfill \\ \hfill & \quad -\frac{2\mu }{3}\left(\frac{1+\nu }{1-\nu }\right)\omega (r).\hfill \end{align} \tag{ 74 }$$

Due to the symmetry of the problem, we have σ_ϕϕ (r) = σ_θθ (r), as can be readily confirmed by a direct calculation. Similarly to the radial component of elastic stress, both hoop components of stress vanish identically in the limit where ω(r) = const or, in other words, where defects are distributed spatially homogeneously throughout the volume of the shell.

The hydrostatic pressure developing in the shell as a result of accumulation of defects is

$$\begin{equation}p(r)=-\frac{1}{3}{\sigma }_{ii}(r)=-\frac{1}{3}\left[{\sigma }_{rr}(r)+{\sigma }_{\phi \phi }(r)+{\sigma }_{\theta \theta }(r)\right].\end{equation} \tag{ 75 }$$

Substituting expressions (73) and (74) in this equation, we find a surprisingly simple relation

$$\begin{equation}p(r)=\frac{4\mu }{9}\left(\frac{1+\nu }{1-\nu }\right)\left[\omega (r)-\bar{\omega }\right].\end{equation} \tag{ 76 }$$

Equation (76) shows that pressure is positive where the local density of relaxation volumes is higher than its mean value, and negative where the local density of relaxation volumes is lower than the mean. The mean pressure, integrated over the entire volume of the shell, is zero.

For the purpose of illustrating the effect of different spatial distributions of defects in the volume of a component, as well as comparing the analytical and FEM solutions, we consider two functions ω(r) plotted in figure 5(a), where case 1 corresponds to equation (C3) and case 2 corresponds to equation (C4) in the appendix. The shell has the inner radius of R₁ = 3.0 m, the outer radius of R₂ = 3.3 m and is assumed to be made of ferritic steel (μ = 80 GPa, ν = 0.29). This is a simplified representation of the vacuum vessel of a fusion reactor. The FEM mesh used in simulations involves 1.33 × 10⁶ hexahedral linear elements. No external stress is applied to the shell.

High energy (14 MeV) fusion neutrons attenuate in steels over distances of the order of tens of cm [22, 23]. For example, in experiments [71] only 10% of the dose from a beam of 14 MeV neutrons reached 40 cm into steel and just 1% of that transmitted dose was comprised of fast (⩾1 MeV in energy) neutrons. Case 2 is representative of these conditions and assumes that ω(r) decays as r⁻² and that a proportion of the neutrons fully penetrate through the 30 cm shell. Case 1, on the other hand, represents a more extreme example of neutrons moving through a highly attenuating material, causing so much damage that swelling saturates at shallow depths, and then being completely stopped within a distance of ∼20 cm. Radiation swelling reaches a dynamic saturation with the density of relaxation volumes of 2% near the plasma, and drops to zero over the distance of approximately 20 cm. The two profiles are chosen so that they have the same total relaxation volume of defects Ω_tot.

The two distributions of ω(r) give rise to very different elastic stress patterns, illustrated in figure 5(b). In both cases, the elastic hoop stresses are of the order of 1 GPa and are higher than the radial component of elastic stress. The hoop stresses are negative close to the inner surface and positive near the outer surface of the shell. Intuitively this is clear since the part of the shell more exposed to irradiation attempts to expand, but is constrained by the less irradiated material, which in turn is under tension to achieve mechanical equilibrium. The maximum value of the von Mises stress is reached at the inner surface and is of the order of 1015 MPa in case 1 and 745 MPa in case 2. Radial displacement fields are also different due to the difference in ω(r), as shown in figure 5(c). However, u_r (r) takes exactly the same values at the inner and at the outer surfaces of the shell. This agrees with equation (71) since Ω_tot is the same in case 1 as in case 2.

A comparison between figures 5(b) and (c) illustrates an important point: a structure that upon a superficial external examination exhibits exactly the same swelling, at least in terms of how its external surfaces move, can develop very different internal stresses. In the examples investigated here, the maximum hoop stress found in case 2 is less than half of that found in case 1.

Contour plots of the internal hydrostatic pressure are compared in figure 6. In both cases, pressure is high and positive close to the inner surface of the shell, vanishes near the centre of the shell, and becomes negative towards the outer surface. The latter may have important implications for the stress-driven diffusion of interstitial elements such as hydrogen, helium or carbon, which would deplete the inner region and segregate towards the outer region of the component. Taking case 1 as an example, the equilibrium concentration of C at 300 K would be about 30% higher on the outside and about 25% lower on the inside than in the unstressed condition.

Zoom In Zoom Out Reset image size

In this section we consider an example where the geometry of a component is too complex to admit an analytical solution, but which is not too dissimilar to what might be considered in the context of a design of a fusion power plant. We analyse, using FEM, the stress fields in a module of a fusion breeding blanket, assuming that it is exposed to a spatially varying flux of neutrons.

The breeding blanket is one of the most significant nuclear components of a fusion power plant, providing and enabling power extraction, tritium fuel sustainability, and radiation shielding [72, 73]. A breeding blanket consists of individual modules containing tritium breeding materials and provides channels for coolant circulation. As an example, we consider a cubic structure with linear dimensions of 500 mm, subdivided into nine submodules. The wall thickness is assumed to be 15 mm, and all the internal edges are filleted with a radius of 3 mm. The structural material for the breeding blanket module is ferritic–martensitic steel Eurofer [74]. The structure is exposed to neutron irradiation from one of the sides, with the direction of the neutron flux being normal to the direction of flow of coolant. In this geometry, radiation exposure varies as a function of only one spatial coordinate [23]. We assume that the maximum value of ω(x) is 2% at the plasma-facing side of the module, see figure 7 of reference [19], and that it decreases down to 0.5% across the component. The geometry of the module and the distribution of ω(x), represented by equation (C5), are shown in figure 7.

Zoom In Zoom Out Reset image size

The component is free of external loads, and is weakly constrained at the four corners of the face z = 0. This is done to prevent its rigid motion but at the same time allows for its free deformation. 2.01 × 10⁶ quadratic tetrahedral elements were used in numerical calculations, following the analysis of convergence given in appendix D.

For the given combination of component geometry and irradiation profile shown in figure 8, there are two locations where stresses are particularly high. Close to the points where internal walls intersect, the local von Mises stress reaches ∼695 MPa. For comparison, the yield point of Eurofer is 530 MPa [75]. At the junction between internal and external walls on the plasma-facing side, the maximum principal stress is close to ∼390 MPa, suggesting that under irradiation it is this location in the breeding blanket module structure that would be susceptible to the nucleation of cracks.

Zoom In Zoom Out Reset image size

From this analysis, and the results given in the preceding sections of this study, we find that the effect of irradiation-induced swelling should be taken into account already at the stage of component design. To exemplify this, we start from the initial design shown in figure 7 and apply two adjustments illustrated in figure 9. First, the internal walls can be rotated by 45° while keeping their length constant. This reduces the number of internal junctions near the plasma-facing side of the component from two to one. Second, the radius of the fillets can be increased where the stresses due to irradiation are expected to be higher. To explore the effect of structural changes, the radius of fillets in the most stressed internal junction was increased to 20 mm, whereas the radius of the four fillets near the plasma-facing side of the component was increased to 30 mm. This produces a modified design, which we subsequently tested using FEM numerical analysis. The von Mises and maximum principal stresses developing in the second design, under the same irradiation conditions, are shown in figure 9, to be compared with figure 8. The FEM model involved 2.11 × 10⁶ elements to ensure convergence. High stresses arise at locations that are similar to those of the first design, but the local von Mises and maximum principal stresses are now lower by 25% and 40%, respectively. The highest maximum principal stress has now developed at other locations in the structure. Further modifications of the design can now be explored in the framework of iterative numerical design optimisation, which is beyond the scope of this work.

Zoom In Zoom Out Reset image size

We observe that irradiation gives rise to sizeable deformations of the component structure. Given the distribution of relaxation volumes of defects, the module expands by approximately 2.5 mm in the x, 3.3 mm in the y and 3.1 mm in the z direction. This expansion is the same for the two design variants considered above, and is close to the deformation of the same structure not containing any inner walls under the same irradiation conditions. The modified module expands by 3.1, 3.4 and 3.3 mm in the x, y and z directions, respectively.

What we find is different from what is known in conventional structural engineering where deformation of a structure subjected to an external load can be reduced by means of internal reinforcements. This methodology does not apply here, since the reinforcements themselves expand following their exposure to irradiation.

Our final note in this section refers to the distribution of hydrostatic pressure p(x) in an irradiated component. Equation (76) shows that in a spherical shell able to expand freely, where both surfaces are subject to the traction-free boundary conditions, p(x) vanishes exactly after averaging over the entire volume of the component. The same applies to the case of a breeding blanket module, provided that it is free of external constraints. The volume-average p(x), evaluated numerically using all the integration points of the FEM model, is of the order of 1 MPa for the design shown in figure 8 and 2 MPa for the modified design shown in figure 9. These values are very small in comparison with the maximum and minimum values of pressure of 285 MPa and −172 MPa computed for the initial component design, and 208 MPa and −139 MPa computed for the modified design, respectively.

The above analysis shows that in the breeding blanket as well as in some other cases that we studied, the yield point of the material was locally exceeded for an assumed swelling profile. To model the resulting response of the structure, the evolution of an irradiated material has to be modelled using elasto-plastic constitutive rules in addition to the laws describing the generation of defects by irradiation. In this context, we expect that the irradiation-driven stress relaxation [76, 77] is going to occur even at relatively lower stresses, partially alleviating them but resulting in the even greater deformation of irradiated components and thus potentially generating stresses in other parts of the reactor structure. Modelling the direct stress-induced plasticity and the slower irradiation-driven stress relaxation, as well as taking into account the temperature and stress dependence of ω_ij (x) defines a broad scope for future work, extending beyond this study. We note that, irrespectively of the occurrence of stress relaxation, the relations between the relaxation volume density ω_ij (x) and elastic stress and strain established in this study remain unchanged, since stress relaxation only affects the form of ω_ij (x) and not the elastic field that it generates.

In this study, we developed a method for computing elastic stress and strain fields developing in structural components of a fusion power plant exposed to irradiation (or indeed any other nuclear facility where components are exposed to a significant irradiation fluence). The approach is based on a defect eigenstrain theorem (equation (8)) that states that the spatially varying density of relaxation volume tensors of defects produced by neutron irradiation (equation (5)), which is a quantity that can be computed from microscopic considerations [30, 31, 33, 47, 48], is identical to the spatially varying field of eigenstrain [39].

In the absence of external geometric constraints or, equivalently, under the traction-free boundary conditions, the elastic stress and internal pressure generated by irradiation always vanish after integration over the entire volume of a component. As a result, a spatially varying exposure of a component to irradiation invariably produces a spatially varying field of stress, compressive at locations where the density of defects is higher, and tensile at locations where the density of defects is lower, see e.g. equation (76). The occurrence of irradiation embrittlement [78] makes operating critical parts of a component under compressive stress preferable to operating under the tensile, negative pressure, conditions.

The magnitude of elastic stress developing in a component as a result of irradiation can be fairly high. It can be estimated as $\sigma (\mathbf{x})\sim \mu [\bar{\omega }-\omega (\mathbf{x})]$, where μ is the shear modulus of the material and $\bar{\omega }-\omega (\mathbf{x})$ is a measure of deviation of the spatially varying density of relaxation volumes of defects from its average value. The density of relaxation volumes of defects is a dimensionless quantity and, depending on the experimental conditions, it varies from a fraction of a percent [25] to several percent [19, 20]. Bearing in mind that for structural steels and tungsten μ is close to a 100 GPa, the irradiation-induced stress can reach several 100 MPa, as illustrated in figure 9. At the same time, a spatially homogeneous distribution of defects does not generate any stress in a free body, although it could still produce large deformations. High stresses can still develop if deformations are constrained.

There are now firm indications suggesting that stress developing in a component as a result of irradiation can elastically polarise defects, resulting in a highly non-linear response of microstructure to stress and irradiation. While this does not alter the fundamental link, equations (9)–(13), between the eigenstrain of the defects and the elastic stress and strain fields that they produce, the combined non-linear self-action of stress and high temperature on the microstructure, already noted in recent simulations and observations [25, 40, 41, 57], requires extensive analysis.

Concluding this study, we point out that the link between the microscopic properties of defects and the macroscopic stress and strain fields that they produce in reactor components, provides an example of a complete multiscale study involving a full treatment of microscopic and macroscopic scales and enabling the assessment of performance of materials and components under realistic operating conditions that can now be quantitatively modelled using advanced supercomputer facilities.

We thank I. Chapman for the provision of laboratory facilities, and are grateful to R. Akers, P.-W. Ma, A. Davis, A. London, F. Hofmann, G. Pintsuk and M. Rieth for helpful discussions. SLD would like to thank A.P. Sutton for a detailed discussion of scientific analysis by J.D. Eshelby [53], which stimulated this study. 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 and 2019–2020 under Grant Agreement No. 633053 and from the RCUK Energy Programme, Grant No. EP/T012250/1. To obtain further information on the data and models underlying the paper please contact PublicationsManager@ukaea.uk. The views and opinions expressed herein do not necessarily reflect those of the European Commission. We gratefully acknowledge the provision of computing resources by the IRIS (STFC) Consortium.

In the following derivation, subscripts after a comma denote differentiation (f_,l = ∂f/∂x_l ), with primed subscripts referring to differentiation with respect to the primed variable (${f}_{ ,{l}^{\prime }}=\partial f/\partial {x}_{l}^{\prime }$).

Consider an elastic body with the volume defined by region V with a surface S = ∂V free of tractions. Mura [39] defines the total strain ${{\epsilon}}_{ij}^{(\text{tot})}$ as the sum of elastic strain _ij and eigenstrain ${{\epsilon}}_{ij}^{\ast }$:

$$\begin{equation}{{\epsilon}}_{ij}^{(\text{tot})}={{\epsilon}}_{ij}+{{\epsilon}}_{ij}^{\ast },\end{equation} \tag{ A1 }$$

with the total strain being compatible ${{\epsilon}}_{ij}^{(\text{tot})}=\frac{1}{2}({u}_{i,j}+{u}_{j,i})$. Similarly, we may introduce the eigenstress ${\sigma }_{ij}^{\ast }$ as a quantity related to eigenstrain ${{\epsilon}}_{ij}^{\ast }$ through Hooke's law:

$$\begin{equation}{\sigma }_{ij}^{\ast }={C}_{ijkl}{{\epsilon}}_{kl}^{\ast }.\end{equation} \tag{ A2 }$$

The Cauchy stress, or elastic stress, is then given by

$$\begin{equation}{\sigma }_{ij}={C}_{ijkl}{{\epsilon}}_{kl}={C}_{ijkl}\left({{\epsilon}}_{kl}^{(\text{tot})}-{{\epsilon}}_{kl}^{\ast }\right),\end{equation} \tag{ A3 }$$

from which follows the equilibrium condition inside the body x ∈ V:

$$\begin{equation}{\sigma }_{ij,j}={C}_{ijkl}\left({{\epsilon}}_{kl,j}^{(\text{tot})}-{{\epsilon}}_{kl,j}^{\ast }\right)=0.\end{equation} \tag{ A4 }$$

The above expression yields a relation between the body force and the eigenstrain:

$$\begin{equation}{C}_{ijkl}{{\epsilon}}_{kl,j}^{(\text{tot})}={C}_{ijkl}{{\epsilon}}_{kl,j}^{\ast }={\sigma }_{ij,j}^{\ast }\equiv -{f}_{ i}.\end{equation} \tag{ A5 }$$

Equivalently, for a body free of external surface forces, the equilibrium condition at the surface x ∈ S follows as

$$\begin{equation}{C}_{ijkl}{{\epsilon}}_{kl}^{(\text{tot})}{n}_{j}={C}_{ijkl}{{\epsilon}}_{kl}^{\ast }{n}_{j}={\sigma }_{ij}^{\ast }{n}_{j}\equiv {t}_{i},\end{equation} \tag{ A6 }$$

where n_j (x) is the outwards pointing normal vector of surface S at point x ∈ S.

The displacement field in a body subjected to a body-force is given by [reference [45], equation (4.12)]

$$\begin{align}\hfill {u}_{j}(\mathbf{\text{x}})& = {\int }_{V}{G}_{jk}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){f}_{ k}({\mathbf{\text{x}}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime }\hfill \\ \hfill & \quad +{\int }_{S}{G}_{jk}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){t}_{k}({\mathbf{\text{x}}}^{\prime }){\mathrm{d}}^{2}{x}^{\prime }\hfill \\ \hfill & \quad -{\int }_{S}{G}_{jk,{p}^{\prime }}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{u}_{i}({\mathbf{\text{x}}}^{\prime }){n}_{m}\enspace {\mathrm{d}}^{2}{x}^{\prime },\hfill \end{align} \tag{ A7 }$$

where f_k is the body force and t_k = C_kpmi u_i,m' n_p is the surface traction.

Now consider a case where the only acting body and surface forces are caused by eigenstrain. We insert the definitions of body force (A5) and traction (A6) into (A7), arriving at:

$$\begin{align}\hfill {u}_{j}(\mathbf{\text{x}})& =-{\int }_{V}{G}_{jk}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{{\epsilon}}_{mi,{p}^{\prime }}^{\ast }({\mathbf{\text{x}}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime }\hfill \\ \hfill & \quad +{\int }_{S}{G}_{jk}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{{\epsilon}}_{mi}^{\ast }({\mathbf{\text{x}}}^{\prime }){n}_{p}\enspace {\mathrm{d}}^{2}{x}^{\prime }\hfill \\ \hfill & \quad -{\int }_{S}{G}_{jk,{p}^{\prime }}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{u}_{i}({\mathbf{\text{x}}}^{\prime }){n}_{m}\enspace {\mathrm{d}}^{2}{x}^{\prime }.\hfill \end{align} \tag{ A8 }$$

We first apply the product rule for differentiation to the first term in (A8), noting that G_jk,p' = −G_jk,p ,

$$\begin{align}\hfill & -{\int }_{V}{G}_{jk}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{{\epsilon}}_{mi,{p}^{\prime }}^{\ast }({\mathbf{\text{x}}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime }\hfill \\ \hfill & \qquad =-{\int }_{V}{\left[{G}_{jk}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{{\epsilon}}_{mi}^{\ast }({\mathbf{\text{x}}}^{\prime })\right]}_{,{p}^{\prime }}{\mathrm{d}}^{3}{x}^{\prime }\hfill \\ \hfill & \qquad \quad -{\int }_{V}{G}_{jk,p}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{{\epsilon}}_{mi}^{\ast }({\mathbf{\text{x}}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime },\hfill \end{align} \tag{ A9 }$$

and then proceed with applying the divergence theorem, to obtain

$$\begin{align}\hfill & -{\int }_{V}{G}_{jk}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{{\epsilon}}_{mi,{p}^{\prime }}^{\ast }({\mathbf{\text{x}}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime }\hfill \\ \hfill & \qquad =-{\int }_{S}{G}_{jk}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{{\epsilon}}_{mi}^{\ast }({\mathbf{\text{x}}}^{\prime }){n}_{p}\enspace {\mathrm{d}}^{2}{x}^{\prime }\hfill \\ \hfill & \qquad \quad -{\int }_{V}{G}_{jk,p}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{{\epsilon}}_{mi}^{\ast }({\mathbf{\text{x}}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime }.\hfill \end{align} \tag{ A10 }$$

Substituting (A10) into (A8), we arrive at

$$\begin{align}\hfill {u}_{j}(\mathbf{\text{x}})& =-{\int }_{V}{G}_{jk,p}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{{\epsilon}}_{mi}^{\ast }({\mathbf{\text{x}}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime }\hfill \\ \hfill & \quad +{\int }_{S}{G}_{jk,p}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{u}_{i}({\mathbf{\text{x}}}^{\prime }){n}_{m}\enspace {\mathrm{d}}^{2}{x}^{\prime }.\hfill \end{align} \tag{ A11 }$$

In preparation of the next step, it is helpful to express the second term in (A11) as a volume integral by using the divergence theorem, noting that G_jk,pm' = −G_jk,pm,

$$\begin{align}\hfill & {\int }_{S}{G}_{jk,p}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{u}_{i}({\mathbf{\text{x}}}^{\prime }){n}_{m}\enspace {\mathrm{d}}^{2}{x}^{\prime }\hfill \\ \hfill & \qquad ={\int }_{V}{\left[{G}_{jk,p}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{u}_{i}({\mathbf{\text{x}}}^{\prime })\right]}_{,{m}^{\prime }}{\mathrm{d}}^{3}{x}^{\prime }\hfill \\ \hfill & \qquad =-{\int }_{V}{G}_{jk,\mathrm{p}\mathrm{m}}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{u}_{i}({\mathbf{\text{x}}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime }\hfill \\ \hfill & \qquad \quad \times {\int }_{V}{G}_{jk,p}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{u}_{i,{m}^{\prime }}({\mathbf{\text{x}}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime }.\hfill \end{align} \tag{ A12 }$$

Substituting (A12) into (A11), we find

$$\begin{align}\hfill {u}_{j}(\mathbf{\text{x}})& =-{\int }_{V}{G}_{jk,p}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{{\epsilon}}_{mi}^{\ast }({\mathbf{\text{x}}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime }\hfill \\ \hfill & \quad -{\int }_{V}{G}_{jk,\mathrm{p}\mathrm{m}}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{u}_{i}({\mathbf{\text{x}}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime }\hfill \\ \hfill & \quad +{\int }_{V}{G}_{jk,p}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }){C}_{kpmi}{u}_{i,{m}^{\prime }}({\mathbf{\text{x}}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime }.\hfill \end{align} \tag{ A13 }$$

Defining the volume-averaged distortion tensor ${\bar{u}}_{j,l}$ by the relation

$$\begin{equation}V{\bar{u}}_{j,l}={\int }_{V}\enspace {u}_{j,l}(\mathbf{\text{x}}){\mathrm{d}}^{3}x,\end{equation} \tag{ A14 }$$

we substitute (A13) into (A14) and write the volume integral over d³ x' as a convolution

$$\begin{align}\hfill V{\bar{u}}_{j,l}& =-{\int }_{V}{C}_{kpmi}\left({G}_{jk,pl}\ast {{\epsilon}}_{mi}^{\ast }\right)(\mathbf{\text{x}}){\mathrm{d}}^{3}x\hfill \\ \hfill & \quad -{\int }_{V}{C}_{kpmi}\left({G}_{jk,pml}\ast {u}_{i}\right)(\mathbf{\text{x}}){\mathrm{d}}^{3}x\hfill \\ \hfill & \quad +{\int }_{V}{C}_{kpmi}\left({G}_{jk,pl}\ast {u}_{i,{m}^{\prime }}\right)(\mathbf{\text{x}}){\mathrm{d}}^{3}x.\hfill \end{align} \tag{ A15 }$$

We proceed with simplifying the second line in (A15). From the definition of the elastic Green's function, we know that [Sutton [45], equation (4.8)]

$$\begin{equation}{C}_{kpim}{G}_{ij,mp}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime })=-{\delta }_{jk}\delta (\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }),\end{equation} \tag{ A16 }$$

which, after manipulation of indices and making use of the symmetries of the stiffness tensor and the Green's function, is equal to

$$\begin{equation}{C}_{kpmi}{G}_{jk,mp}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime })=-{\delta }_{ji}\delta (\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime }).\end{equation} \tag{ A17 }$$

The substitution of (A17) into the second line of (A15) yields

$$\begin{align}\hfill & -{\int }_{V}{C}_{kpmi}\left({G}_{jk,pml}\ast {u}_{i}\right)(\mathbf{\text{x}}){\mathrm{d}}^{3}x\hfill \\ \hfill & \qquad =-{\int }_{V}{\int }_{V}\left[-{\delta }_{ji}{\delta }_{,l}(\mathbf{\text{x}}-{\mathbf{\text{x}}}^{\prime })\right]{u}_{i}({\mathbf{\text{x}}}^{\prime }){\mathrm{d}}^{3}{x}^{\prime }\enspace {\mathrm{d}}^{3}x\hfill \\ \hfill & \qquad ={\int }_{V}{u}_{j,l}(\mathbf{\text{x}}){\mathrm{d}}^{3}x=V{\bar{u}}_{j,l},\hfill \end{align} \tag{ A18 }$$

where we used ∂_x (δ*f)(x) = ∂_x f(x) = f'(x). Cancelling out $V{\bar{u}}_{j,l}$ in (A15) hence leads to

$$\begin{align}\hfill & {\int }_{V}{C}_{kpmi}\left({G}_{jk,pl}\ast {{\epsilon}}_{mi}^{\ast }\right)(\mathbf{\text{x}}){\mathrm{d}}^{3}x\hfill \\ \hfill & \qquad ={\int }_{V}{C}_{kpmi}\left({G}_{jk,pl}\ast {u}_{i,{m}^{\prime }}\right)(\mathbf{\text{x}}){\mathrm{d}}^{3}x.\hfill \end{align} \tag{ A19 }$$

Next, substituting

$$\begin{equation}{C}_{kpmi}{u}_{i,{m}^{\prime }}={C}_{kpmi}{{\epsilon}}_{im}^{(\text{tot})}={C}_{kpmi}({{\epsilon}}_{im}+{{\epsilon}}_{im}^{\ast }),\end{equation} \tag{ A20 }$$

into the second line of (A19) and simplifying the resulting expression, we arrive at

$$\begin{equation}0={\int }_{V}{C}_{kpmi}\left({G}_{jk,pl}\ast {{\epsilon}}_{im}\right)(\mathbf{\text{x}}){\mathrm{d}}^{3}x.\end{equation} \tag{ A21 }$$

The above expression is in effect a volume integral over a convolution between two functions f and g, which may also be expressed as

$$\begin{equation}{\int }_{V}\left(f\ast g\right)(\mathbf{\text{x}})\enspace {\mathrm{d}}^{3}x=\left({\int }_{V}f(\mathbf{\text{x}}){\mathrm{d}}^{3}x\right)\left({\int }_{V}g(\mathbf{\text{x}}){\mathrm{d}}^{3}x\right).\end{equation} \tag{ A22 }$$

Applying equation (A22) to equation (A21), we arrive at

$$\begin{equation}0={C}_{kpmi}\left[{\int }_{V}{G}_{jk,pl}(\mathbf{\text{x}}){\mathrm{d}}^{3}x\right]V{\bar{{\epsilon}}}_{im},\end{equation} \tag{ A23 }$$

where ${\bar{e}}_{im}={V}^{-1}\int {{\epsilon}}_{im}(\mathbf{\text{x}}){\mathrm{d}}^{3}x$ is the volume-averaged elastic strain tensor.

The so-called auxiliary tensor

$$\begin{equation}{D}_{jkpl}={\int }_{V}{G}_{jk,pl}(\mathbf{\text{x}}){\mathrm{d}}^{3}x,\end{equation} \tag{ A24 }$$

is not generally zero, enabling us to conclude that the volume-averaged elastic strain tensor vanishes

$$\begin{equation}{\bar{{\epsilon}}}_{im}={V}^{-1}{\int }_{V}{{\epsilon}}_{im}(\mathbf{\text{x}}){\mathrm{d}}^{3}x=0.\end{equation} \tag{ A25 }$$

Note that the same applies to the volume-averaged elastic stress tensor, ${\bar{\sigma }}_{im}=0$.

Finally, recalling the definition of the volume relaxation tensor

$$\begin{equation}{{\Omega}}_{ij}={\int }_{V}{{\epsilon}}_{ij}^{(\text{tot})}(\mathbf{\text{x}}){\mathrm{d}}^{3}x=V{\bar{{\epsilon}}}_{ij}^{(\text{tot})}=V{\bar{{\epsilon}}}_{ij}+V{\bar{{\epsilon}}}_{ij}^{\ast },\end{equation} \tag{ A26 }$$

and substituting ${\bar{{\epsilon}}}_{ij}=0$, we arrive at the central result

$$\begin{equation}{{\Omega}}_{ij}=V{\bar{{\epsilon}}}_{ij}^{\ast }.\end{equation} \tag{ A27 }$$

To conclude, the elastic field _ij does not contribute to volume change in a body with the surface free of tractions. Any volume change is solely effected by the volume-average eigenstrain field ${{\epsilon}}_{ij}^{\ast }$, as shown by expression (A27). The proof above is general insofar as the eigenstrain field may be non-zero across the whole of the body, including the surface.

We can also recover the well known relation between relaxation volume tensor and dipole tensor [8, 44]

$$\begin{equation}{{\Omega}}_{ij}=V{\bar{{\epsilon}}}_{ij}^{\ast }={S}_{ijkl}{C}_{klmn}V{\bar{{\epsilon}}}_{mn}^{\ast }={S}_{ijkl}{P}_{kl},\end{equation} \tag{ A28 }$$

where ${P}_{kl}=V{\bar{\sigma }}_{kl}^{\ast }$ is the dipole tensor.

Consider a body that underwent deformation such that point x in the initial configuration is at x' in the deformed configuration. The initial and deformed coordinates are related through the displacement field

$$\begin{equation}{\mathbf{\text{x}}}^{\prime }=\mathbf{\text{x}}+\mathbf{\text{u}}(\mathbf{\text{x}}).\end{equation} \tag{ B1 }$$

Considering the displacement as a coordinate transformation from x → x', the volume of the body is given by the volume integral

$$\begin{equation}V={\int }_{V}\mathrm{det}(\mathbf{J})\mathrm{d}x\enspace \mathrm{d}y\enspace \mathrm{d}z,\end{equation} \tag{ B2 }$$

where det(J) is the determinant of the Jacobian matrix

$$\begin{equation}{J}_{ij}=\frac{\partial {x}_{i}^{\prime }}{\partial {x}_{j}}={\delta }_{ij}+{u}_{i,j}(\mathbf{\text{x}}),\end{equation} \tag{ B3 }$$

given by

$$\begin{align}\hfill \mathrm{det}(\mathbf{J})& =1+{u}_{i,i}+\frac{1}{2}\left({u}_{i,i}{u}_{j,j}-{u}_{i,j}{u}_{j,i}\right)\hfill \\ \hfill & \quad +{u}_{1,i}{u}_{2,j}{u}_{3,k}{{\epsilon}}_{ijk},\hfill \end{align} \tag{ B4 }$$

where _ijk is the Levi-Civita tensor.

In the limit of small strain

$$\begin{equation}\mathrm{det}(\mathbf{J})\approx 1+{u}_{i,i},\end{equation} \tag{ B5 }$$

leading to the expression for the volume of the body [49]

$$\begin{equation}V={\int }_{V}\left(1+{u}_{i,i}\right)\mathrm{d}x\enspace \mathrm{d}y\enspace \mathrm{d}z.\end{equation} \tag{ B6 }$$

With this in mind, even if the displacement field has been obtained using linear elasticity theory, expression

$$\begin{equation}{{\Omega}}_{\mathrm{t}\mathrm{o}\mathrm{t}}={\int }_{V}{{\epsilon}}_{ii}^{(\text{tot})}(\mathbf{x})\mathrm{d}V,\end{equation} \tag{ B7 }$$

only captures the volume change of the distorted body to first order in strain. This is typically a valid approximation as mean strains in linear elasticity theory are usually limited to a few percent.

In the treatment of distortions in a rectangular foil considered in section 4, function ω(z) was taken to be proportional to the ion implantation profile. Three tungsten ion implantation profiles were considered in section 4, corresponding to ion energies of 10, 20 and 50 MeV. The data points that were generated using the SRIM package [58, 59] are shown in figure 1. These three data sets were fitted using the function

$$\begin{equation}f(\zeta )=\frac{{c}_{1}+{c}_{2}\zeta +{c}_{3}{\zeta }^{2}}{{\left[1+\mathrm{exp}\left(\frac{\zeta -{c}_{4}}{{c}_{5}}\right)\right]}^{{c}_{6}}},\end{equation} \tag{ C1 }$$

where ζ = h − z is the depth variable used in figure 1. Using the least squares regression we determined the fitting parameters given in table 1.

Function ω(ζ) was then assumed to be proportional to f(ζ). The proportionality constant was set to 0.034 37 so that the peak of the 20 MeV curve corresponds to the relaxation volume density of 1%.

For the cylindrical tube case explored in section 5, the profile of ω(ρ) plotted in figure 4(a) is

$$\begin{equation}\omega (\rho )=\frac{0.01}{{(4.0-\rho )}^{2}},\end{equation} \tag{ C2 }$$

where 2.5 < ρ < 3.0 is in mm.

For the case of a spherical shell described in section 6, the profiles of ω(r) plotted in figure 5(a) are

$$\begin{equation}\omega (r)=\frac{0.02}{1+\mathrm{exp}\left(\frac{r-3.15}{0.01}\right)},\end{equation} \tag{ C3 }$$

for case 1 and

$$\begin{equation}\omega (r)=\frac{0.005}{{(r-2.5)}^{2}}-0.002\enspace 7843,\end{equation} \tag{ C4 }$$

for case 2. Both functions are defined on the interval 3.0 < r < 3.3, where r is given in metre units. The second term in equation (C4) ensures that integral ∫_V ω(x)dV has the same value in both cases.

For the breeding blanket module of section 7, function ω(x) is defined as

$$\begin{equation}\omega (x)=\frac{1069}{{(x-731)}^{2}},\end{equation} \tag{ C5 }$$

where 0 < x < 500 is expressed in mm units.

As opposed to the analysis given in the sections preceding section 7, there is no analytical solution for the breeding blanket module described in section 7 to assess the quality of the FEM mesh. To obtain the numerical FEM results presented in section 7, we started from a relatively coarse mesh of about 5 × 10⁵ tetrahedral elements and applied to the entire model the adaptive remeshing algorithm of Abaqus, reaching a final mesh of about 2 × 10⁶ elements. In such a way, the programme was allowed to decrease the element size down to about 0.1 mm at locations where the gradients of stress were relatively high, at the same time coarsening the mesh where gradients were relatively small.

Figure 10(a) shows the highest values of the von Mises and maximum principal stresses as a function of the total number of elements included in the FEM representation of the first design with horizontal and vertical inner walls. The mesh refinement procedure was stopped after two iterations as the change in these highest stresses was considered to be acceptably small (3% for the von Mises, 1% for the maximum principal stress). A similar mesh refinement process was carried out for the modified design model with slanted inner walls and increased fillet radii, as shown in figure 10(b). Since the region where the maximum principal stress was evaluated (i.e. the inset of figure 9(b)) was no longer the maximum for the entire structure, the second iteration caused the value to decrease slightly. The numerical values of stress converged to within 1% both for the von Mises and for the maximum principal stress.

Zoom In Zoom Out Reset image size