Finite element models for radiation effects in nuclear fusion applications

Deuterium-tritium fusion reactions produce energy in the form of 14.1 MeV neutrons, and hence fusion reactor components will be exposed to high energy neutron irradiation while also being subjected to thermal, mechanical and magnetic loads. Exposure to neutron irradiation has numerous consequences, including swelling and dimensional changes, comparable in magnitude to the peak transient thermal deformations occurring in plasma-facing components. Irradiation also dynamically alters the various thermo-mechanical properties, relating temperature, stress and swelling in a strongly non-linear way. Experimental data on the effect of neutron exposure spanning the design parameter space are very sparse and this highlights the relevance of computer simulations. In this study we explore the equivalence between the body force/surface traction approach and the eigenstrain formalism for treating anisotropic irradiation-induced swelling. We find that both commercial and massively parallelised open source software for finite element method (FEM) simulations are suitable for assessing the effect of neutron exposure on the mechanically loaded reactor components. We demonstrate how two primary effects of irradiation, radiation swelling and the degradation of thermal conductivity, affect the distributions of stress and temperature in the divertor of the ITER tokamak. Significant uncertainties characterising the magnitude of swelling and models for treating it suggest that on the basis of the presently available data, only an order of magnitude estimate can be given to the stress developing in reactor components most exposed to irradiation during service.


Introduction
Several countries have set an ambitious goal of delivering energy from fusion power plants to the grid by the 2040s.To meet this timescale, exploiting the capabilities of virtual engineering is going to be a necessity [1] and the strategy Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
for the DEMO reactor needs to move beyond the design-byexperiment approach [2].Burning deuterium-tritium plasma is going to be producing energy in the form of 14.1 MeV neutrons, and yet performing neutron irradiation experiments at the component scale with a fusion neutron energy spectrum, reaching the high dose expected in an operating reactor, is still not possible.Experimental facilities presently under construction either ignore neutron irradiation effects [3] or focus on testing materials solely on a small spatial scale [4].
The low availability of irradiation facilities capable of delivering the relevant high dose, the difficulty of extracting in-situ microstructural information spanning the reactor design-relevant range of stress and temperature operating conditions, and the dose rate effect implying the need to conduct experiments on the timescale comparable with components' lifetime, all suggest the desirability of developing robust mathematical and computational models to estimate the operating lifetime and explore the failure scenarios of reactor components ahead of its construction.Furthermore, exploring a virtual fusion reactor model offers the most natural way of defining the operating parameter space for materials, including temperature, irradiation exposure and stress, for the subsequent investigation in a dedicated materials testing facility [4,5].
Properties of materials operating in a radiation environment change as a function of dose, dose rate, temperature, the initial composition and microstructure of a material, and other factors.Neutrons alter the chemical composition of materials through transmutations, as has been extensively documented for tungsten [6,7].Nonelastic neutron-nuclei reactions give rise to the intense generation of γ-photons in mediumto-heavy elements [8,9], which contribute to components' internal volumetric heating in addition to the direct heating by the plasma.
Materials change dimensions and shape due to the combined effect of swelling and creep [10,11].Swelling spans the range from a fraction of a percent to several percent in steels and structural alloys and is markedly sensitive to temperature [10][11][12][13][14]. Swelling and creep can be anisotropic due to the anisotropy of the crystal structure [15][16][17] or due to the evolution of radiation defects that relax the applied stress [12,[18][19][20].Moreover, lattice defects and transmutation products elastically distort the lattice and act as scattering centres for electrons and phonons [21], thereby degrading the thermal conductivity [22].
Structural design criteria for the in-vessel components of ITER mandate the consideration of irradiation effects [23].While the neutron exposure of structural components in ITER is expected to be lower than in DEMO or in a commercial power plant, the total accumulated dose is expected to be many times the end of life dose in the pressure vessel of a fission power plant [24], making it necessary to treat radiation exposure in the design of a tokamak.For example, volumetric swelling is a significant source of elastic stress developing in reactor components exposed to neutrons [25,26].
Since performing irradiation tests on the components scale under fusion neutron energy spectra is still not feasible, quantifying effects of irradiation in silico can be pursued as an alternative [27].A large amount of experimental data on neutron irradiated materials are available from fission reactors or spallation neutron sources.These are expensive tests, lasting over several months, which are still not able to explore the high-dose limit and the specific effects of the 14 MeV neutron energy spectrum characterising fusion irradiation.Computer simulations, on the other hand, suffer from the lack of quantitative representation of fundamental microstructural phenomena and hence still involve relatively large uncertainties.
An extensively studied component of fusion reactors such as ITER or DEMO is a 'monoblock' tile consisting of a CuCrZr inner pipe, cooling the plasma-facing W armour, with a Cu interlayer joining the cooling pipe with the armour.Although several finite element method (FEM) studies of a monoblock were performed, involving coupled thermomechanical analysis [28][29][30][31][32], less attention was devoted to the investigation of effects of irradiation [2], even though the divertor is expected to be exposed to neutron fluxes that are among the highest anywhere in the DEMO reactor structure [8].The central issue with the structural analysis of a monoblock is the lack of detailed materials data.This, as well as the recognition of uncertainties inherent in the models and the relatively limited predictive strength of simulations in the treatment of fracture phenomena [33] are among the points highlighted below.
Structural simulations of fusion components involve the treatment of neutron and γ-radiation transfer and interaction with materials, magneto-thermo-mechanical loads, and plasma and fluid dynamics [34].Interaction of neutrons with materials gives rise to phenomena spanning a broad range of length-and time-scales.Recognizing this fact, FEM developed for nuclear fission applications is multiscale and computationally massively parallelized [35,36].In fusion, high performance computing (HPC) has been deployed in computationally intensive materials simulations [37], in plasma physics turbulence simulations [38], as well as in neutron transport simulations [39].It has not yet been as widely used in structural mechanics.In the context of FEM simulations, HPC enables extending the method in two directions.First, it increases the number of elements accessible to the treatment, translating the new quantitative capabilities into greater geometric complexity and higher fidelity of the simulations.Second, it enables performing high-throughput simulations [40,41] to quantify uncertainties of the results, given the inherent scatter of input parameters or material properties affected by irradiation.
Bearing in mind the objective of implementing fusion neutron radiation effects in HPC-FEM, we use the MOOSE (multiscale object-oriented simulation environment) software [42], already extensively used in the context of nuclear fission applications [36,43].In this study we do not focus on large-scale calculations, but rather concentrate on the methodology and comparison with the simulations performed using the commercial FEM code Abaqus, demonstrating the consistency of the two methods.
In what follows, we further develop and apply the mathematical methodology for treating radiation swelling in FEM.We compare approaches where it is treated either as eigenstrain [25,44], or as a distribution of body forces and surface tractions [27].Subsequently, we explore case studies pertinent to fusion engineering: a single crystalline material exhibiting anisotropic swelling under applied external stress, a spherical shell as a greatly simplified model for a vacuum vessel, an anisotropic ion-irradiated zirconium sample and an ITER divertor monoblock.For some of the examples, exact analytical solutions are given and compared with numerical results.

The linear elastic problem in the presence of radiation defects
Strain is a quantity that measures the deformation of a material.There are two kinds of strain: elastic strain and eigenstrain.The former is deformation that is linked to the presence of elastic stress.The latter is deformation of other physical origin, and since it can exist in the absence of an applied stress, it is also called a stress-free strain.Common examples of eigenstrain are thermal expansion, phase transformations and plastic deformation [44,45].Irradiation-induced swelling is another source of eigenstrain that is of particular importance to nuclear applications [25].Below, we consider thermal expansion and radiation swelling.
Although radiation swelling is often assumed to be a scalar quantity, it is in general anisotropic [12].Imagine a general configuration of defects produced by irradiation, where each defect has a tensorial relaxation volume Ω (a) kl [27].Its trace Ω (a) kk characterises the amount of volumetric expansion or contraction caused by the defect.For instance, Ω (a) kl for a vacancy is isotropic and the relaxation volume of a vacancy is negative, but for a dislocation loop it is anisotropic, with positive or negative relaxation volume depending on its interstitial or vacancy nature [46].To quantify the collective effect of an ensemble of defects, we define the density of relaxation volume tensors [25,27] where each defect is located at R a , δ() is the Dirac delta function and x = (x, y, z) ≡ (x 1 , x 2 , x 3 ) is the spatial coordinate.Despite being conceptually defined in terms of isolated defects, the density of relaxation volume tensors can be generalised to dislocation loops and networks of dislocations [47].
As shown in [25], the dimensionless tensorial swelling ω kl (x), treated as a continuous differentiable function, has the same properties as eigenstrain.In particular, its trace characterizes the degree of local expansion or contraction of the material if the material is unconstrained.Let us start from the total strain that can be directly measured using spatially-resolved x-ray or electron diffraction [48,49], taking in this case the name of lattice strain.This is the sum of elastic strain ε kl and any source of eigenstrain ε * kl , in the small strain limit where this additive decomposition remains valid [45].The total strain is found by the differentiation of the displacement vector field u(x) [50], where an index after a comma indicates partial differentiation with respect to a component of coordinate vector x with that index.In what follows, we take eigenstrain in the form where α kl is the thermal expansion tensor and T is the temperature relative to a reference temperature state.Elastic stress σ ij is computed by combining equation ( 2) and the tensor of elastic constants c ijkl as where the three tensorial quantities c ijkl , α kl (x) and ω kl (x) in general can also depend on temperature [51].Symmetrising the derivative of displacements in ( 5) is no longer necessary as the tensor of elastic constants is symmetric with respect to the transposition of its last two indices [50].If the eigenstrain distribution is known, calculating the resulting stress amounts to solving a direct eigenstrain problem [45], and this is what we focus on below.
The material can be subjected to body forces due to gravity or magnetism.Both loads are present in a breeding blanket of a fusion reactor, which is close to strong magnets and where the candidate materials are ferromagnetic ferritic steels such as Eurofer97.In a static analysis, Cauchy's equation of mechanical equilibrium reads where ρ is the mass density, g i the acceleration of gravity and f m i the magnetic body force [52,53].
If we assume that parameters c ijkl and α kl are spatially homogeneous, we arrive at a linear differential equation for the vector field of displacements u(x): This equation also shows that since the sources of eigenstrain do not contain the displacement variable, the effect of spatially varying eigenstrain is the same as that of effective body forces [25,44].
If the elastic behaviour, thermal expansion, and swelling are isotropic, i.e.
where µ is the shear modulus, ν is the Poisson ratio, δ ij is the Kronecker delta and α is the linear thermal expansion coefficient, equation ( 7) takes a simpler form, where we drop the explicit spatial dependence to simplify the notation, A solution can be obtained for a given set of boundary conditions.Typically, these are of mixed type.For example, vector u(x) can be directly prescribed for the constrained surfaces.Other surfaces may have a known distribution of externally applied stress or temperature.In the absence of external loads, the surfaces are traction free, i.e. σ ij (x)n j (x) = 0 at all the surfaces, where n(x) is the external surface normal vector at point x.In the presence of eigenstrains, applying traction free conditions to equation ( 5) implies that at surfaces Eigenstrains that are non-zero at the surface of a body act as effective surface tractions [44].
The total force acting on a component of volume V, as a result of body forces that also capture the effect of swelling, is The total force acting on its external surface ∂V, as a result of the tractions that capture the effect of the same swelling, is Recalling the divergence theorem for a generic tensor we see that the effective body and surface forces always form a mechanically self-equilibrating system.

Equivalence of the eigenstrain and body force/surface traction formalisms
We now focus purely on radiation swelling, demonstrating how it can be numerically treated in FEM in two different but entirely equivalent ways, namely as an effective thermalexpansion-like eigenstrain or as a combination of body forces and surface tractions.
In equation ( 4) thermal expansion and swelling linearly add up in the definition of eigenstrain.It is possible to make use of this fact and represent swelling by a fictitious thermal expansion.In the most general form, this requires a spatially dependent tensor of thermal expansion coefficients.If one is not solving a coupled thermo-mechanical problem, this is the most straightforward option.However, it is not convenient if temperature is also part of the analysis and if the material properties are temperature-dependent, although using the thermalexpansion-like analogy is still possible [2].In this case it is advisable to use the equivalence between an arbitrary eigenstrain and a suitable distribution of body forces f i and surface tractions t i .Recalling equations (2.12) and (2.13) by Mura [44], these distributions are which in the elastically isotropic case become This can be further simplified if swelling is isotropic, namely where B = 2µ 3 If swelling is represented by a fictitious thermal expansion, an FEM code treats it as eigenstrain.In other words, FEM evaluates elastic strain by subtracting eigenstrain from the total (lattice) strain, see equation ( 2).
One disadvantage of using the distribution of body forces (15) and tractions ( 16) is that these are not forces of true mechanical origin like gravity, but effective forces representing the effect of eigenstrain.Both contribute to the total displacement field in the body and to the stress field that we call σ tot ij .However, this total stress, unlike elastic stress, is not an observable but rather a formal quantity given by the tensor product of elastic constants and the total strain.We must subsequently subtract the eigenstress contribution from this total stress, just like in equation ( 5), to find the elastic stress The second term in the right-hand side of this equation simply reads Bω(x)δ ij if both tensors are isotropic.
The need for this step can be intuitively understood using the simplest possible example: a free cubic sample of a material with sides L expanding due to homogeneous swelling ω or, equivalently, a homogeneous temperature increase.It is clear that the cube expands homogeneously and then permanently remains under zero stress.The gradient term in equation ( 19) vanishes everywhere, implying that no body forces are present.Still, there are surface tractions following from equation (20), which are of magnitude Bω at all the faces of the cube.This self-equilibrating system of forces produces the desired observed expansion of the cube.However, the tractions also generate a state of spatially homogeneous hydrostatic tension with σ tot ij = Bωδ ij that is constant throughout the internal volume of the cube.The elastic-stress-free condition σ ij (x) = 0 everywhere in the cube is then enforced by subtracting from the computed total stress the quantity Bωδ ij , which is the second term in equation (21).A notable shortcoming of the above body force approach is that the setup can be cumbersome to implement since every surface can have a different spatially-varying distribution of surface tractions.We now provide an example where the numerical values produced using the two approaches are compared.Let us consider a cube that in the presence of defects undergoes anisotropic dimensional changes caused by swelling that spatially varies in magnitude but retains a constant tensorial character, i.e.
For instance, this can occur in the case where irradiation produces [100]-type dumbbell self-interstitial atom defects in a bcc material, and the material is under a uniaxial stress of sufficient magnitude to orient most of these defects along the same direction, in this example with the axes of all the defects pointing along the y coordinate axis.Using the relaxation volume tensor of a [010] dumbbell interstitial in Fe [54], we take that where β ii = 1.The density of relaxation volumes in this example is assumed to have the form where d(x, y, z) is the distance from the centre of the cube, where we show an explicit application of equation (21).We emphasize that the stress evaluated by FEM where the body forces and surface tractions method is employed, is incorrect unless equation ( 21) is subsequently used.

A spherical shell
In this section, we compare the fields of stress and deformations predicted by Abaqus and MOOSE for a case that is also tractable analytically: a spherical shell where the tensorial density of relaxation volumes of radiation defects is directionally isotropic and is only a function of the radial coordinate.
The shell has the internal radius R 1 and external radius R 2 .Thermal deformations are not considered.An exact analytical solution was derived in [25,27] and we simply quote it here.Given the spherical coordinate system (r, θ, ϕ), and assuming that the radial displacement is The three components of elastic stress that can be computed from this radial displacement are )] Here, ω is the volume-average density of relaxation volumes of defects in the material where ) is the volume of the shell.A finite element model for a simplified fusion reactor vacuum vessel was created using one eighth of a spherical shell with R 1 = 7.0 m and R 2 = 7.3 m, made of ferritic steel (E = 210 GPa and ν = 0.3).The density of relaxation volumes, chosen in an analytical form ω(r) = 10 −4 /(r − 6.9) 2 , where r is given in metres, was selected in such a way so that ω(r) would decay by a factor of 10 from its maximum value ω(R 1 ) = 1% at the inner surface of the shell over the distance of ∼20 cm, to mimic the attenuation of the neutron flux in steel.To mesh the geometry, 10 6 hexahedral linear elements were used.Symmetric boundary conditions were applied to the three orthogonal sections.
A quantitative comparison between analytical and numerical solutions is given in figure 2. The three data sets are in excellent agreement.In particular, after averaging over three different linear paths to eliminate the possible mesh effects, the mean difference between the FEM and analytical solutions was 0.6 MPa (von Mises stress) and 0.4 MPa (hydrostatic stress) where MOOSE was used, and 1.1 MPa (von Mises stress) and 0.6 MPa (hydrostatic stress) in the case of Abaqus calculations.As the magnitude of the stress components themselves is of order of hundreds of MPa, the relative error of the hydrostatic stress is similar for the two codes, being 0.8% for MOOSE and 1.0% for Abaqus.
For the chosen profile of ω(r), equation ( 26) predicts the outward radial displacements of 5.727 mm and 5.973 mm for the internal and external surfaces, respectively.The corresponding values calculated by Abaqus are 5.731 mm and 5.977 mm, whereas MOOSE gives 5.730 mm and 5.976 mm.The numerical values are very close and the relative error with respect to the exact analytical result equation ( 26) is less than 0.1%.shows ω(x) in the cross-section at z = 100 mm and β ij is given by equation (23).ω ij can be described either by using fictitious thermal expansion, i.e. the setup shown in (b) producing results in (d), or by a combination of body forces and tractions (BF/T), with the setup in (c) giving results in (e).The displacement field is the same for both setups, whereas if the thermal expansion setup produces the elastic stress σ ij , BF/T setup yields the total stress σ tot ij .If the thermal expansion setup is used, swelling is already represented as a (thermal) eigenstrain.If the body forces and tractions method is used, we must remember that the elastic stress is evaluated as σ ij = σ tot ij − c ijkl ω kl , see also equation (21).To illustrate that the two stress tensors differ by the amount c ijkl ω kl , the von Mises stress and the hydrostatic stress are plotted along different linear paths in (f ) and (g) respectively.The mesh of the one eighth of the shell that was modelled using the three symmetry planes is shown below the legend, with the indication of the location of the three radial paths on which the stress was extracted.Although slight differences due to mesh effects are visible, they are very small compared with the magnitude of the stresses.

Ion-irradiated anisotropic zirconium foil
In this section we explore effects of swelling in a crystallographically anisotropic material, using ion-irradiated Zr an example.The crystal lattice of Zr has the hexagonal close packed structure, characterised by a considerable degree of elastic anisotropy.Moreover, the swelling of Zr following its exposure to irradiation is also highly anisotropic [15,16,55].To account for the swelling anisotropy, we express the eigenstrain in a form where β ij is assumed to have distinct eigenvalues.This form of the eigenstrain holds if the exposure is spatially dependent but the crystal orientation and the character of swelling remains invariant, like in the case of a single-crystal thin film.
In the low temperature limit, experiments show that the Zr lattice expands isotropically and two-dimensionally in the basal plane, whereas it contracts in the c-direction by approximately half of the basal plane expansion [55].On the other hand, simulations found the same trend, but with a smaller degree of contraction in the c-direction, compared with the observed expansion in the basal plane [16].
A finite element model for a 80 µm × 80 µm × 10 µm foil was created using 250 000 elements.The normal to the foil and the c-direction were parallel to the z-axis.One fourth of the foil was modelled exploiting symmetries.The foil was constrained at its bottom surface (u z = 0) and swelling was considered to depend only on z. ω(z) was assumed to have saturated at 1% in the top 2 µm of the sample, decaying to zero with the maximum penetration depth of about 4 µm.The values of elastic constants of Zr are c 11 = c 22 = 155 GPa, c 33 = 172 GPa, c 12 = 67 GPa, c 13 = 65 GPa, c 44 = 36 GPa [56].
Molecular dynamics simulations of low-temperature selfion irradiated tungsten, in agreement with observations, found that very early in the irradiation history, the self-interstitial and vacancy defects occur primarily in isolation.Defects of the former type strongly interact [46] and rapidly evolve into a dislocation network [11,47] that grows to the point of forming sheets of interstitials that, as a result of the effect of boundary conditions, are oriented parallel to the foil surface [20].
From these observations, we postulate two forms for the tensorial part of eigenstrain (30).These are meant to model highly idealised scenarios and are the simplified representations of observed swelling.The first case represents a homogeneous dispersion of isolated vacancies and self-interstitials.The second assumes that all the self-interstitials have formed perfect crystallographic planes orthogonal to the z direction, which is therefore the only direction in which they expand the material, with additionally the same homogeneous dispersion of vacancies.Based on ab initio values for the relaxation volumes of self-interstitials and vacancies in Zr [16], in the 'isolated defect' case we assume On the other hand, if complete crystallographic planes are formed, we obtain The justification for the tensorial forms above is given in appendix A. These forms were chosen intentionally so that the traces of the two tensors are identical and positive.In the 'isolated defect' case, the irradiated material attempts to expand against the resistance from the unirradiated layer underneath.In the case where extra complete crystallographic planes form through the coalescence of self-interstitial defects, the opposite effect occurs, but the magnitude of the in-plane mismatch is lower.Hence, we expect that the irradiated layer should be under hydrostatic compression in the former and hydrostatic tension in the latter case, with the opposite stress state developing in the unirradiated layer at the bottom of the sample due to mechanical self-equilibration.Stresses are expected to be lower in the case where self-interstitial defects have self-assembled in the form of complete crystallographic planes.
The problem can be readily solved using a suitable FEM code, in our case either Abaqus or MOOSE, producing effectively the same result, as illustrated in figure 3. The simulated elastic stress and displacement fields are consistent with expectations; the maximum absolute magnitude of the hydrostatic stress is 298 MPa in the first and 72 MPa in the second case.The elastic strain energy in the case of isolated interstitials is 3.29 × 10 −4 J mm −3 , and it is 0.19 × 10 −4 J mm −3 in the case of perfectly formed extra atomic planes.In the former case the energy is 17.3 times larger (which incidentally is very close to 17.1, which is the square of the ratio of maximum hydrostatic stresses computed for the two cases).This suggests that there is an elastic driving force for irradiated films that arranges the defects so that the swelling is anisotropically concentrated in the out-of-plane direction.This agrees with experimental observations on ion-irradiated W [20] and SiC [19].

Radiation swelling of an ITER divertor monoblock
To remove the heat deposited by the exhaust plasma at the divertor, which is going to experience the highest heat loads over the entire tokamak reactor structure, the design chosen for ITER involves an array of Cu-alloy pipes shielded by tungsten armour tiles, because of the high melting point and low sputtering erosion rate of the latter [32].The pipes, as per ITER specifications, are made of a precipitation-hardened Cu alloy containing 0.6%-0.9%Cr and 0.07%-0.15%Zr [57], usually referred to as CuCrZr.A soft pure Cu interlayer joins the pipe to the surrounding W material.In addition to the intense heat load approaching 10-20 MW m −2 [31], the divertor plasma-facing components also have to withstand intense neutron irradiation.We therefore include the irradiation effects in the model.
We start by exploring a simplified case that admits an analytical solution: a pipe without the W tile, shown in figure 4(a).In an ITER monoblock, a CuCrZr pipe with thickness of 1.5 mm is joined to a 28 × 28 × 12 mm 3 W tile through a 1 mm thick Cu interlayer.Following [2], we adopt a uniform isotropic swelling model for W (ω W ) and Cu (ω Cu ), which assumes random orientations of radiation defects [27] and neglects the swelling of CuCrZr.The numerical values of swelling ω W = 0.6% and ω Cu = 1.0%have the same ratio as in [2] but are of lower magnitude.Since the W part dominates the volume of the component, and since W material is about four times as stiff as the pipe material, the deformation of the monoblock is dominated by the properties of W alone.
Consider two co-axial fully bonded cylinders under plane strain conditions.The inner surface of the CuCrZr pipe is at R 1 = 6 mm; the contact surface between CuCrZr and the Cu layer is at R 2 = 7.5 mm and the contact surface between Cu and W is at R 3 = 8.5 mm.The radial displacement of the Cu layer and of the CuCrZr pipe, and their common strain along the axis of the pipe are These expressions can now be inserted into equation (5) where, neglecting the effects of thermal expansion in the isotropic elastic case, we find Each of the two materials develops a different stress state.
Denoting µ 1 and ν 1 for Cu and µ 2 and ν 2 for CuCrZr, we find that the non-vanishing stress components are The displacement at the outer boundary of the Cu layer is defined by the swelling of the adjacent W part of the monoblock, and so is the axial displacement, which is the plane strain condition (35).Additionally, the inner surface of the CuCrZr is considered to be traction free.The coolant pressure of 4 MPa would add a negligible correction to the result below and is therefore neglected.We can define the four boundary conditions to solve for the constants a 1 , a 2 , b 1 and b 2 .The boundary conditions and the algebraic expressions for the four constants are given in appendix B.
This simplified cylindrical pipe-only model enables the validation of mechanical contact condition in the FEM framework, but may also be used for a rapid assessment of the effect of variation of swelling in W and in Cu, different geometries and materials.The incorporation of thermal expansion, radiation swelling in CuCrZr and effects of internal pressure is conceptually straightforward, but it was not included here to not encumber the expressions.The analytical solution does not apply if symmetry is broken by the presence of heat flux at the W top surface; but it is still of interest as an exact result describing the case where thermal loads are off.
The mesh of the pipe-only model is shown in figure 4(a), and the results are plotted in figure 4(b) alongside the analytical expressions ( 33), ( 34) and ( 37)- (42).One quarter of the pipe was modelled exploiting symmetries, and very good agreement between analytical and numerical solutions is evident.
We now proceed to simulating a full FEM model for an ITER monoblock.The geometry follows [32]; one quarter of the tile was discretized into 510 000 linear hexahedral elements; the mesh is shown in figure 5.The contact surfaces between CuCrZr and Cu, and between Cu and W are assumed to be fully bonded.The mechanical boundary conditions are the symmetry conditions along the x and z symmetry planes and u y = 0 at the bottom W surface.The thermal boundary conditions are as follows.Heat removal by the water coolant at temperature of 100 • C was simulated by convective heat transfer with a transfer coefficient of 110 kW (m 2 •K) −1 [32].The reference temperature for the evaluation of deformation from thermal expansion was set to 100 • C. Radiative heat dissipation was allowed at external surfaces.Numerical values of materials' properties are given in table 1.The monoblock was exposed to the thermal flux of 20 MW m −2 over a period of 10 s.This is twice the design steady state value in ITER and is expected to be reached during transient load oscillations [31].
At present, experimental tests under irradiation at the component level are unattainable.Although levels of irradiation in ITER are going to be relatively modest, yet many times the end-of-life dose of 0.045 dpa in a fission reactor vacuum vessel [24], the methodology for treating neutron irradiation effects in plasma-facing component design is important for DEMO and for commercial reactors [2,33].With the aim of applying our formalism, in what follows we consider an artificially high degree of exposure.Three key irradiation effects that were included in the analysis are: 1. radiation-induced swelling, 2. transmutation-and defect-induced degradation of thermal conductivity, 3. heating by γ-radiation, as detailed in the following paragraphs.
1.The dimensions of a monoblock tile are small compared to the neutron attenuation distance, which is of the order of 10 cm.Hence, the radiation exposure of the monoblock is approximately constant.A strong swelling gradient due to differential exposure is not expected, but there can still be one due to the swelling versus temperature variation effect.Swelling in tungsten [13], copper and its alloys [14], stainless steels [62] and chromium [63] has a characteristic temperature dependence, showing a maximum of swelling within a certain temperature interval specific to each element.The swelling versus temperature profile is broadly symmetric around the peak swelling temperature, with a maximum value of the order of 1%-10%.For simplicity, in this study we approximate the temperature dependence of swelling by a Gaussian profile For W, we assumed a peak swelling value of 1.0% at 750 • C with τ = 100 • C; for Cu and CuCrZr we assumed a peak swelling value of 2.0% and 0.5%, respectively, at 325 The swelling profile given by equation ( 43) depends only on temperature.On the other hand, in equation ( 1) we introduce a spatially dependent swelling.This highlights an important point, namely that irradiation induces stresses either if a gradient of exposure or of temperature or, of course, of both is present.In the specific example of the monoblock, exposure to irradiation is broadly constant; it is the enormous temperature gradient that can introduce sharp gradients in swelling.
Here, we have limited the analysis to the case of isotropic swelling [27] and excluded any stress relaxation or plasticity phenomena.These are beyond the scope of this study, and require a careful ad hoc treatment and hence cannot be appropriately described by conventional thermal creep or plasticity based on gliding dislocations.However, both are important [29,31] and should be included when extending this initial study.
2. The accumulation of radiation defects and transmutation products causes a notable decrease in thermal conductivity, especially in W, as it was recognised in earlier studies [2,64].This affects the temperature profile in the monoblock and, as a consequence, the stress from thermal expansion as well as swelling.Thermal annealing of the defects leads to an increase in the thermal conductivity of irradiated tungsten as temperature increases, contrary to the decrease in this parameter in an unirradiated material.
3. Neutron irradiation is a source of volumetric heating mediated by the colossal generation of γ-photons in the medium to heavy elements like W and Cu.Considering the DEMO neutron flux at the divertor region, we estimated the volumetric heating power of about 3 Watts per gram, on the basis of the analysis carried out in [8].
We performed two transient heating simulations, one considering the scenario not involving irradiation and one introducing the three radiation-induced effects: swelling, reduced thermal conductivity, and the gamma heating.Similarly to the approach taken by Fursdon and You [2], we included the effect of irradiation on thermal conductivity and added radiation swelling to the setup of the simulation; we additionally also included the constant volumetric thermal power due to the neutron-induced gamma emission [8,9].The results corresponding to the final time step are shown in figure 6.
Table 1.Material properties used in the ITER monoblock simulations described in section 6.The parameters vary as functions of temperature T: thermal conductivity κ, specific heat cp, Young's modulus E, Poisson's ratio ν, density ρ, the linear thermal expansion coefficient α.We note that both radiation swelling and the degradation of thermal conductivity are expected to arise over a time scale that is many orders of magnitude longer than the 10 s pulse; the pulse is therefore thought to occur after these irradiated conditions have already been attained.The neutron exposure producing swelling of the order of a percent [13] is greater than that foreseen after the entire duration of ITER operations.The irradiation scenario is not a novel failure pathway for the ITER divertor, but a demonstration of the simulation methodology and an approximation for the effects of long-term exposure of a monoblock to fusion neutrons, as foreseen in a steady-state DEMO operation mode.Although the exact geometry of the monoblock will likely be different, the materials selected for the DEMO divertor are currently the same [2,33].
The above irradiation scenario is as a representation of the steady-state DEMO operation.We note that the number of pulses necessary for accruing the dose of several dpa would probably cause failure due to thermo-mechanical fatigue well before the onset of substantial swelling.We used a simple mapping from the temperature, recorded at every nodal point, to local swelling and changes in thermal conductivity, simplifying the history dependence of microstructure on variable temperature, which warrants a separate analysis.
The simulations taking into account irradiation effects show substantially higher stress levels due to the additional source of eigenstrain associated with radiation defects, as well as higher temperatures because of the lower thermal conductivity of irradiated materials.The maximum temperature at the corners of the top surface increased from 2005 • C to 2305 • C. The average temperature of the top surface increased from 1870 • C to 2160 • C. The volume average temperature of the upper half of the W material increased from 860 • C to 1025 • C, or by about 20%.The temperature of the CuCrZr pipe was affected to a lesser extent; the average temperature of the upper half volume, for instance, increased by 5% from 360 • C to 380 • C. Of the two main effects of irradiation on the temperature profile, namely the volumetric gamma heating and the reduction in thermal conductivity, the latter has a substantially greater effect.This was confirmed by repeating the simulation without the body heat source and finding the temperatures only a few percent lower.
The stress profile was completely altered by the presence of swelling.In W, radiation swelling is most apparent in the temperature range from 600 • C to 900 • C [13].The highest temperature at the W/Cu interface was about 640 • C, rising continuously towards the top surface, and hence leading to the formation of a spatial interval between the pipe and the top surface where the maximum of swelling was attained across the width and the depth of the tile.This maximum swelling band is very visible in the von Mises map shown in figure 6(b).Under the unirradiated conditions, the maximum principal stress in W was the highest, at 1470 MPa, at the interface with Cu at the centre of the tile and at the side closest to the top surface.Under irradiation, it rose to 1835 MPa at a different location along the W/Cu interface, close to where the dashed line in figure 6(a) is.Similarly, the von Mises stress peaked at the top of the W/Cu interface at the outer surface in the unirradiated case, but with swelling included, the maximum was higher by about 20% and moved close to the diagonal direction.The CuCrZr/Cu pipe was even more affected by irradiation.The maximum principal stress in the unirradiated case was the highest in CuCrZr, attaining the value of 130 MPa, with stress in Cu being close to 100 MPa.With radiation swelling included, the stress level of 485 MPa was reached in Cu, with stress in CuCrZr rising to 460 MPa.Such a stress state would clearly generate substantial plastic deformation especially in the pipe materials.In this work, we did not consider the effect of irradiation on plastic behaviour and fracture toughness, which are the two serious concerns in relation to structural integrity [33].With respect to previous work we confirmed the role of Values of the maximum principal stress and temperature are also plotted in (c) along a path, indicated by a dashed line in (a), running at 3 mm depth parallel to the front surface.The CuCrZr/Cu and Cu/W interfaces are at 1.5 mm and 2.5 mm respectively.In (a), thermal expansion is the only source of stress.In (b), swelling is additionally present and is responsible for the high stress where the temperature conditions favour the maximum of radiation swelling.Moreover, temperature is higher in (b) as a result of the irradiation-degraded thermal conductivity.Hollow symbols are the Abaqus data and filled symbols are the MOOSE data; colour scales in (a) and (b) are identical for comparison.variation of thermal conductivity [2], and additionally focused on spatially-varying and temperature-dependent swelling.To address lifetime predictions in the future, in components that must endure higher exposure than ITER, changes in plastic and thermal and swelling behaviours ought to be considered.
As a final note, exploring the trends rather than providing exact numerical values was the goal of this section; predictive calculations require going beyond the linear elastic limit and the isotropic swelling approximation.What the calculations do show, however, is that experimentally verified values of swelling can completely change the results.As long as the uncertainties associated with swelling and irradiation creep are large, both in terms of experimental observations and the underlying models, the design of plasma facing reactor components cannot be fully validated and certified.The predicted stress levels strongly depend the peak swelling parameter ω max chosen in our simulations to be close to 1.0%, 2.0% and 0.5% for W, Cu and CuCrZr in equation (43).These values vary, depending on chemical composition of the material, and hence are sensitive to transmutations, to the radiation exposure as well as to the neutron energy spectra, the dose rate and other factors.Hence, the calculated in-service stress levels can presently be regarded only as estimates until better data and models become available to represent the stresssensitive anisotropic swelling in materials under the expected fusion irradiation conditions.
To further illustrate the point, we performed a brief survey of the swelling parameter space.Experiments [13,14] often find higher values of swelling than what we used so far in simulations.Accordingly, in previous FEM models the maximum swelling values were set to 3.2% and 5.0% in W and Cu respectively [2].We then ran four more simulations setting ω max in W in the range between 0.5% and 2.0%, ω max in Cu between 0.5% and 4.0%, and keeping everything else unchanged, including the constant value of 0.5% for ω max in CuCrZr.The effect of these parameters on the stress state in the monoblock was significant and the results are summarised in table 2 using the highest value of the maximum principal stress found in the three materials as the assessment criterion.
The results show that this particular indicator referring to the pipe materials is particularly sensitive to swelling in Cu.This appears significant as the Cu interlayer was selected with the intention of reducing the interfacial stress.The highest maximum principal stress in W appeared less clearly linked to swelling in either W or Cu, but generally it increased with increasing swelling.Not only the magnitude but also the location of the highest maximum principal stress value changed.In W it was located close to the Cu/W interface, but for instance if ω max was chosen at 2.0% in both materials, the highest value was found along the external edge of the tile, where the peak swelling temperature was reached.

Conclusions
Anisotropic irradiation-induced swelling can be included in FEM either by using a distribution of body forces and surface tractions, or by directly treating swelling as a spatially varying field of eigenstrain.The open-source code MOOSE and the commercial code Abaqus were employed to study the stress distribution in a neutron-irradiated spherical shell, in an ion-irradiated zirconium and in an ITER monoblock.
During the plasma heating exposure of the monoblock under the heat load of 20 MW m −2 , the maximum thermal expansion gradient of approximately 1% over 6 mm distance develops between the plasma-facing surface and the cooling pipe.Once the steady-state conditions of neutron and thermal exposure are reached during the foreseen DEMO operationsalthough admittedly the steady state is unlikely to be reached in ITER-the irradiation can generate a distribution of eigenstrain that is comparable or higher in magnitude than thermal expansion, but over a shorter distance of about 1 mm, if driven by comparable temperature gradients.
Uncertainties in the swelling input data and constitutive models dominated all the other sources of uncertainty for the monoblock simulations; reducing the former is therefore very important for a reliable engineering design.For instance, our FEM simulations showed that radiation swelling in the Cu interlayer joining tungsten to CuCrZr could make the presence of this interlayer a detrimental, rather than beneficial, factor in the design due to the interfacial stress build-up between the W plasma facing part of the monoblock and the cooling pipe.One of the ways to address the uncertainty in the data is running parametric analyses.Completing the ten time steps of the monoblock calculations of section 6 required of the order of 100 core-hours.Repeating similar simulations for hundreds or even thousands of times, as well as increasing the number of elements to faithfully represent larger structures, would rapidly evolve into a major supercomputing application as opposed to a task suitable for a conventional desktop computer.Up-scaling FEM to HPC is therefore necessary in the context of realistic fusion reactor design studies.
Swelling, temperature and stress are closely coupled and the success of testing reactor components in silico depends sensitively on the reliability of an underlying model describing how these three quantities affect one another under reactor operating conditions, as well as on their temporal history coupled with the history of microstructural evolution.

2 ) 2 .
The elastic constants are set to µ = 80 GPa and ν = 0.3.The length of the side of the cube is L = 20 cm and the external stress σ yy = 100 MPa.The two different FEM setups are shown in figures 1(a)-(c).The resulting stress and deformations are given in figures 1(d)-(g),

Figure 1 .
Figure 1.A cubic sample of iron swells anisotropically following exposure to irradiation, see section 3 for further detail: the simulation setup is in (a)-(c) and results are in (d)-(g).Stress and deformation are determined from the spatially varying density of relaxation volumes of defects ω ij (x) = ω(x)β ij ; where (a)shows ω(x) in the cross-section at z = 100 mm and β ij is given by equation(23).ω ij can be described either by using fictitious thermal expansion, i.e. the setup shown in (b) producing results in (d), or by a combination of body forces and tractions (BF/T), with the setup in (c) giving results in (e).The displacement field is the same for both setups, whereas if the thermal expansion setup produces the elastic stress σ ij , BF/T setup yields the total stress σ tot ij .If the thermal expansion setup is used, swelling is already represented as a (thermal) eigenstrain.If the body forces and tractions method is used, we must remember that the elastic stress is evaluated as σ ij = σ tot ij − c ijkl ω kl , see also equation(21).To illustrate that the two stress tensors differ by the amount c ijkl ω kl , the von Mises stress and the hydrostatic stress are plotted along different linear paths in (f ) and (g) respectively.

Figure 2 .
Figure 2. Spherical shell deformed by a spherically symmetric swelling field.The FEM results (discrete symbols) and the analytical solution are in very good agreement, as shown by plotting the von Mises stress (a) and the hydrostatic stress (b) along a radial path from the inner to the outer shell surface.The mesh of the one eighth of the shell that was modelled using the three symmetry planes is shown below the legend, with the indication of the location of the three radial paths on which the stress was extracted.Although slight differences due to mesh effects are visible, they are very small compared with the magnitude of the stresses.

Figure 3 .
Figure 3. Distribution of hydrostatic stress in an anisotropically swollen Zr film shown in (a) and (b) for the cases of isolated self-interstitials and perfect crystallographic planes formed by the coalescence of these defects, respectively.The full agreement found between the two FEM codes is illustrated in (c) as a function of depth from the irradiated top surface.The foil either expands or contracts in the c-direction, confirmed by the plot of uz, as a result of the anisotropy.Hollow symbols are the Abaqus data and filled symbols are the MOOSE data.

Figure 4 .
Figure 4. (a) Model of a simplified monoblock model where it is assumed that the surrounding W material imposes boundary conditions on the CuCrZr/Cu pipe (dimensions in mm); symmetry was enforced on planes normal to x and z.(b) plots of the stress components and radial displacement resulting from the setup in (a).The analytical expressions of section 6 are compared to FEM.Hollow symbols are the Abaqus data and the filled symbols are the MOOSE data.

Figure 5 .
Figure 5. Mesh and geometric parameters of an ITER monoblock, consisting of a CuCrZr inner pipe and a Cu interlayer at the centre of a W tile. Dimensions are in mm.The two symmetry planes that are exploited are normal to the x and to the z directions.

Figure 6 .
Figure 6.The Von Mises stress (shown in the left-hand side) and temperature (shown in the right-hand side) are compared for the (a) un-irradiated and (b) irradiated case at the end of a 10 s interval of transient heating by the exposure of the top surface to the heat flux of 20 MW/m 2 .Values of the maximum principal stress and temperature are also plotted in (c) along a path, indicated by a dashed line in (a), running at 3 mm depth parallel to the front surface.The CuCrZr/Cu and Cu/W interfaces are at 1.5 mm and 2.5 mm respectively.In (a), thermal expansion is the only source of stress.In (b), swelling is additionally present and is responsible for the high stress where the temperature conditions favour the maximum of radiation swelling.Moreover, temperature is higher in (b) as a result of the irradiation-degraded thermal conductivity.Hollow symbols are the Abaqus data and filled symbols are the MOOSE data; colour scales in (a) and (b) are identical for comparison.

Table 2 .
(43)est values of the maximum principal stress σ MP in an ITER monoblock evaluated for various peak swelling values in equation(43)for Cu and W. Swelling in CuCrZr was kept constant at 0.5%.The location of the peak value of the maximum principal stress is also affected by the choice of ωmax.