Multi-scale fracture probability analysis of tungsten monoblocks under fusion conditions

Plasma facing components inside future nuclear fusion reactors are subjected to a high heat load and intense irradiation conditions. Using advanced computational material models, several problems can be solved that reflect tungsten monoblocks under fusion relevant loading scenarios. This allows for the identification of the conditions under which material failure is probable. The material model and parameters are identified such that the mechanical behaviour is in accordance with the homogenized behaviour of a previously developed crystal plasticity model on the microscopic scale. The heterogeneous stress field that follows is analysed in order to assess the probability of material failure, which is typically reflected by unstable crack propagation. Since fracture is an inherently multi-scale problem, critical regions are analysed in detail by means of a representative volume element. The resulting analysis reveals that in case the stress relaxation in the monoblock under the applied static heat load is complete, the probability of unstable crack propagation can reach values close to 35%. Finally, the impact of prolonged neutron irradiation is simulated by means of a cluster dynamics model. Although irradiation drastically increases the brittleness of tungsten, its impact on the overall monoblock performance remains limited.


Introduction
Nuclear fusion has the potential to help in solving the world's energy problem. This is partly due to the abundant availability of its required fuel. In fusion reactors, the hydrogen isotopes deuterium and tritium are fused to form helium. During this process, a small amount of mass is converted into a vast amount of energy. The extremely high temperature that is required for this to occur causes the particles involved to behave as a plasma. The plasma facing components in the reactor are subjected to extreme conditions, including an exceptionally high heat flux and intense neutron irradiation.
In the reactor design, pure tungsten monoblocks are incorporated which should provide the required protection against these extreme conditions. While the highly demanding conditions are known to significantly deteriorate the mechanical properties of tungsten, a considerable lifetime of these components is mandatory for fusion reactors to be economically viable [1]. The monoblocks are subjected to a staggering heat load and simultaneously water cooled. The strong temperature gradient that is hereby induced, results in a spatially dependent thermal expansion and hence a significant stress field.
The mechanical response at the length scale of the monoblock is governed by the underlying microstructure, which consists of an assembly of grains. Within these grains, plastic deformation takes place * Corresponding author. along a limited number of energetically favourable crystallographic slip planes. Assuming a pre-existing distribution of microscopically small cracks or flaws in the tungsten microstructure, Griffith's criterion [2] reveals that unstable crack propagation emanates from a high level of stress. This cleavage type of fracture initiates at a very small length scale, and rapidly propagates to macroscopic dimensions pertinent to that of the monoblock.
In this work, the monoblock is considered to be fully recrystallized. The corresponding texture consists of crystal orientations with a uniformly random distribution. In order to describe the mechanical response at the macroscopic scale, an isotropic Von Mises type of plasticity model is employed that captures the homogenized stress response of a crystal plasticity model developed in earlier work. The current work analyses critical regions of the monoblock in detail by prescribing the deformation and temperature evolution of a macroscopic material point to a representative volume element that resembles the underlying microscopic scale.
In earlier work [3], the influence of fusion-relevant neutron irradiation on the mechanical response of a tungsten microstructure is investigated by means of cluster dynamics simulations [4]. It was shown that tungsten becomes increasingly brittle upon neutron irradiation, as its brittle-to-ductile transition temperature increases. This aspect is also taken into account in the current analysis, where the impact on the performance of a monoblock is evaluated for various levels of neutron irradiation.
The staggering heat flux in combination with the irradiation conditions with a highly specific neutron spectrum [5] are unique to nuclear fusion conditions and are particularly challenging to reproduce experimentally. The numerical models are based on scarce experimental data and are therefore vital in order to assess the performance of the tungsten components under fusion conditions. A novel aspect of the current work is the multi-scale analysis of tungsten monoblocks under both an extreme heat load as well as neutron irradiation. This framework includes an innovative multi-scale study of the probability of unstable crack propagation.
The outline of this paper is as follows. First, in Section 2 the relevant microscopic models are discussed. These include the crystal plasticity material model and cluster dynamics model adopted from earlier work [3]. Then in Section 3, the macroscopic boundary value problem is discussed, along with the material model used at this length scale. In Section 4, the failure probability model is discussed, which applies to both length scales and is based on the weakest-link theory. The results for heating and cooling of the monoblock are presented in Section 5. Finally, the main findings of the current work are discussed and summarized in Section 6.

Microscopic model
In order to describe the mechanical behaviour of polycrystalline tungsten, a representative volume element (RVE) is utilized. It consists of a grid-like geometry of 4 × 4 × 4 grains, discretized into 12 × 12 × 12 elements. The RVE is visualized in Fig. 1 where the colours indicate the orientation in height direction. In the figure, HD, TD and WD resemble the height direction, tube direction and width direction, respectively. The RVE is subjected to 3D periodic boundary conditions. The grains have random crystal orientations, such that the homogenized behaviour of the RVE can be considered isotropic. More details on the real recrystallized tungsten microstructure can be found in literature [6] and also our earlier work [7].
In order to describe the mechanical material response of tungsten at the microscopic scale, a crystal plasticity model is presented in Section 2.1. In this model, the overall plastic deformation consists of crystallographic slip on energetically favourable planes. The resistance to slip is known to increase upon neutron irradiation, as it induces lattice defects that impede the motion of dislocations. The evolution of irradiation induced defects is simulated by means of a cluster dynamics model in Section 2.2.

Crystal plasticity model
In order to describe the material behaviour at a material point on the microscopic scale, a finite strain crystal plasticity model is adopted which is based on the work of Asaro [8]. Since this model allows to take into account the orientation of crystallographic slip systems, the constitutive relation is inherently anisotropic. Furthermore, the resistance to slip is strongly dependent on temperature and defect concentrations. The crystal plasticity material model is implemented in the finite element software Marc Mentat using element type 57 and a user subroutine.

Kinematics
The deformation gradient tensor describes the mapping of a material point from the reference configuration to the deformed configuration. It is decomposed into an elastic part, a thermal part and a plastic part: with e the elastic deformation gradient tensor and p the plastic deformation gradient tensor. Due to this decomposition, an intermediate configuration can be defined in which only the plastic deformation is applied. The thermal expansion is isotropic and is described by means of the thermal deformation gradient tensor th as where is the difference between the initial temperature and the current temperature. Furthermore, is the second-order unit tensor. The plastic velocity gradient tensor is given by wherėis the slip rate on slip system , with = 1, … , s . Vectors ⃗ 0 and ⃗ 0 are the slip direction and slip plane normal of slip system in the undeformed (or intermediate) configuration. The Schmid tensor 0 = ⃗ 0 ⃗ 0 is given by their dyadic product. Following earlier work [9], only the {1 01 } ⟨111⟩ family of slip systems is assumed to be potentially active [10][11][12].

Elastic contribution
The second Piola-Kirchhoff stress tensor in the intermediate configuration is given by the double inner product between the fourth-order stiffness tensor and the elastic strain tensor: where 4 C is the fourth-order isotropic elasticity tensor, defined by the Poisson's ratio and the temperature dependent Young's modulus . The elastic Green-Lagrange strain tensor e in the intermediate configuration is given by where e = T e ⋅ e is the elastic right Cauchy-Green deformation tensor. Note that e can easily be rewritten into other stress measures as well, for example: Here, is the second Piola-Kirchhoff stress tensor in the initial configuration and the Cauchy stress tensor in the current configuration.

Crystallographic slip
The nucleation and propagation of dislocations through the tungsten lattice results in crystallographic slip. Here, a rate-dependent framework is adopted where the slip rate on slip system is phenomenologically described by means of a power-law: The material parameters anḋ0 are the strain-rate sensitivity parameter and reference slip rate, respectively. The resolved shear stress on slip system is given by: The resistance to crystallographic slip is decomposed into three parts: where the obstacle slip resistance obs is further discussed in Section 2.  Here, the grain volume is taken to be 1.64 × 10 6 3 , in line with earlier work [7]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) 2.1.3.1. Obstacle slip resistance. The dislocation density typically increases with plastic deformation, which is associated with an increase in slip resistance. The evolution of the obstacle slip resistance is given by [13]: where ℎ are the hardening coefficients, given by Here, are the coefficients of the interaction matrix. Its diagonal components = 1 resemble the self-hardening coefficients, while for its off-diagonal components = 1.4 is taken. Material parameters ℎ 0 , and ∞ are the initial hardening rate, hardening exponent and saturation slip resistance respectively.

Kink-pair slip resistance.
In order for dislocations to bypass the periodic resistance of the crystal lattice, they have to overcome the thermally activated Peierls barrier. The formation of kink-pairs along the dislocation line is the most energetically favourable way to overcome these barriers. The kink-pair slip resistance is given by [14]: Here, kp,0 is the kink-pair slip resistance at = 0 K. Furthermore, exponents and are material constants. The expression for the knee-temperature k on slip system is taken from previous work [3]: with the ''knee-slip-rate''̇k =̇k p,0 exp(− 0 ∕( B )) [3], B the Boltzmann constant, and 0 the activation energy. Finally, = 10 −7 s −1 is an arbitrary small value that prevents kp to be highly sensitive to low values oḟ< .

Cluster dynamics
During the nuclear fusion process, deuterium and tritium atoms are fused to form a helium atom and a highly energetic neutron. Since neutrons have no charge, they are not magnetically confined in the torus and irradiate the components inside the reactor. A highly energetic neutron can collide with a tungsten atom in one of the plasma facing components. This tungsten atom is then referred to as the primary knock-on atom and displaces many subsequent atoms through the tungsten lattice, leaving a trace of defects. These defects consist of vacancies V and self-interstitials I. Since single defects are mobile, they can interact with one another. Within a very short time frame after a collision, most defects typically annihilate, i.e. V + I → 0. Defects can also grow in size by absorption of a defect of the same type i.e. I +I → I +1 or V +V → V +1 , where denotes the defect cluster size. Frequent occurrence of the latter interaction reaction results in vacancy clusters of a large size, that are also referred to as voids or cavities. In line with the work of Shah et al. [15], the backwards interaction reactions V → V −1 + V and I → I −1 + I are taken into account as well. Defects can also shrink in size according to V + I → V −1 or I + V → I −1 . Finally, defects can also be annihilated due to their interaction with dislocations according to D+ → D, where resembles either vacancies V or self-interstitials I.
All of the above interaction reactions are taken into account by means of the master equations provided in A. Here, the concentrations of vacancy clusters and interstitial clusters are denoted by V and I respectively. The rate at which defects are produced due to neutron irradiation is discussed in Section 2.2.1. The rate at which the interaction reactions occur scale with the coefficients , and . Their superscript + or − indicate whether a forward or backward reaction is considered respectively. These coefficients are dependent on temperature and the defect size considered, as detailed in Section 2.2.2. The formulation of these coefficients is largely adopted from the work of Li et al. [4].

Defect production rate
The rate of the total number of defects, per unit of volume, produced in the tungsten lattice due to neutron irradiation is given by tot . The work of Sand et al. [16] has shown that during defect production, a fraction of the vacancies and interstitials cluster together during a very short time scale. These fractions are indicated by c V and c I , respectively. The rate at which unclustered defects are produced is therefore obtained by the first line of the following equation: for = V, I. Regarding the clustered defects, the work of Sand et al. [16] has shown that their distribution can be accurately described by − . Here, are temperature dependent exponents that describe the distribution of the produced defect clusters. The work of Setyawan et al. [5] also showed that there is a marked upper limit to the size at which defects are produced during neutron irradiation. The expression is therefore adjusted such that no defects are produced of a cluster size c and larger, see Eq. (14). The pre-factors can be computed by equating the second line of Eq. (14) to the product c tot , which is the rate at which clustered defects are produced: The expression is obtained by equating the sum of the defects categorized in each defect size to the total number of clustered defects.

Rate coefficients
In this section, the coefficients that correspond to the interaction reactions between defects are discussed. The formulation of these temperature and defect size dependent coefficients is largely adopted from the work of Li et al. [4].

Rate coefficients for vacancy clusters.
The rate at which the concentration of vacancy clusters of size decreases due to the reaction equation V +V → V +1 is given by the product + V V . Here, the rate coefficient + is given by: where V = (3 ∕4 ) 1∕3 + √ 3 0 ∕4 is the capture radius of a vacancy cluster of size , = 3 0 ∕2 is the atomic volume and 0 is the lattice parameter. The diffusivity of a single vacancy is given by: Here, the pre-exponential factor 0 V and single-vacancy migration energy mig V are material parameters. The above reaction equation can also occur the other way around, i.e. V → V −1 + V. Due to this reaction equation, the concentration of vacancy clusters of size decreases by − V . The rate coefficient − is given by: The binding energy that corresponds to V + V → V +1 is given by Here, b V 2 is the binding energy of a cluster of size 2 and is considered to be a material constant. The formation energy of a single vacancy f V is also considered to be a material constant.

Rate coefficients for interstitial clusters.
Similar to the vacancy clusters discussed in the previous subsection, interstitial clusters can also absorb or eject a single interstitial, i.e. I +I → I +1 and I → I −1 +I. These reactions cause the concentration of interstitial clusters of size to decrease. The corresponding rates are given by + I I and − I respectively. The corresponding coefficients + and − are given by: Similar to Eq. (17), the diffusivity of a single interstitial is described by the Arrhenius equation where material parameters mig I and 0 I are the migration energy and pre-exponential factor. The capture radius of an interstitial loop of size is given by where is the length of the Burgers vector. The capture efficiency I ( = I, V) by an interstitial loop of size is given by [17]: , 1 ) .
The binding energy of an interstitial cluster of size to a single interstitial is given by: where the formation energy of a single interstitial f I is a material parameter, as well as b I 2 which equals the binding energy of an interstitial cluster of size two.

Rate coefficients across species.
The rate at which the concentration of vacancy clusters of size decreases due to the reaction equation V + I → V −1 is given by + V +I V I . Likewise, the rate at which the concentration of interstitial clusters of size decreases due to the reaction equation I + V → I −1 is given by + I +V I V . Coefficients + V +I and + I +V are respectively given by [18]: The reaction equation that corresponds to = 1 (i.e. V + I → 0) is considered to be a special case, and given by Here, IV is a material parameter.

Rate coefficients regarding defect-dislocation interaction.
The rate at which the concentration of a single defect decreases due to the absorption by dislocations (i.e. D + → D) is given by + D+ (with = I, V). The corresponding coefficients + D+ are given by: Here, d is the dislocation density.

Irradiation induced slip resistance
Due to the accumulation of irradiation induced lattice defects, the slip resistance increases. The increase in slip resistance due to a concentration of voids can be described by the dispersed barrier hardening model: is the average spacing between voids of size . Here, the diameter and concentration of a void of size are given by V and V , respectively. Previous work [3] has shown that the influence of neutron irradiation can be adequately described by taking a root-sum-square type of superposition law, which for each slip system gives: Here, is the number defect sizes included in the simulation. Molecular dynamics simulations in literature show that the slip resistance decreases with increasing temperature. Simulations from Terentyev et al. [19] show that this decrease is most pronounced at low temperatures, but saturates at elevated temperatures. Hence, a constant irr is incorporated in order to account for the influence of temperature in the domain that is currently investigated.

Monoblock model
In order to investigate the performance of the monoblocks, their thermo-mechanical response is simulated. The monoblock geometry considered in this work is adopted from Hirai et al. [20] and is shown in Fig. 2, where the dimensions are indicated in mm. The out-of-plane thickness equals 12 mm. In this figure, the grey area resembles tungsten, while orange and yellow represent the CuCrZr cooling tube and the Cu interlayer, respectively. A flow of water is running through the cooling tube causing forced convection.
In Section 3.1, the boundary value problem that corresponds to the analysis of the monoblock is detailed, along with the discretization and boundary conditions used. Then, in Section 3.2, the material model is discussed.

Boundary value problem
In order to reduce the computational cost, two existing planes of symmetry are exploited. This allows for simulating only a quarter of the monoblock geometry, as is shown in Fig. 3. The geometry is discretized into serendipity brick elements with reduced integration, i.e. 20-node elements with 8 integration points. The geometry of the tungsten part is discretized into 1504 elements, while the cooling tube and interlayer consist of 320 and 400 elements, respectively. In the figure, the planes of symmetry are indicated by the blue and green area. The displacements of the nodes on these planes in the direction of the surface normal are suppressed. Since the monoblocks mounted on top of the cooling tube have a spacing of 0.5 mm between them, the cooling tube extends by 0.25 mm in the model geometry. The end of the cooling tube is indicated by the orange area in Fig. 3. In order to enforce symmetry, the axial displacements of the nodes on this surface are tied such that their axial displacement is identical. Finally, a pressure of 3.3 MPa is applied inside the cooling tube, i.e. the purple area in Fig. 3.
The red area resembles the plasma-facing surface, on which a heat load of 10 MW∕m 2 is applied. This causes a rapid temperature increase, starting from a homogeneous stress-free reference state with an initial temperature of 0 = 350 K. The influence of transients is disregarded in the current analysis. On the inside of the cooling tube, forced convection by means of a water flow causes a radial heat flux that is expressed as: Here, wall is the temperature at the inside of the cooling tube and bulk is the temperature of the water. The latter is set equal to the homogeneous stress-free reference temperature, i.e. bulk = 0 = 350 K. Figure 4 in the work of Le et al. [21] shows that a heat transfer coefficient ℎ conv = 5 × 10 4 W∕m 2 K is applicable for a wide range of temperatures up to approximately wall = 520 K, which is not exceeded in the current work. The remaining thermo-mechanical parameters are given in B. In the simulations that follow, thermo-mechanical problems  are solved in which the temperature profiles are obtained at every time step.

Macroscopic material model
The crystal plasticity model in Section 2.1 allows for the description of the mechanical behaviour at the level of the individual crystals. Here, the orientation of the slip systems within the grains is taken into account. At the macroscopic length scale (i.e. the scale of the monoblock), the mechanical response is governed by the homogenized behaviour of the polycrystalline material. However, in this work, the macroscopic behaviour is simplified using a Von Mises plasticity model that shows a highly similar response.
The Von Mises yield criterion follows a similar temperature and strain rate dependency as the slip resistance in the crystal plasticity model. Therefore, the following hardening law is adopted in the current work: which is analogous to Eq. (9). Here, obs is included to describe the influence of the obstacle slip resistance on the yield stress. To this end, the formulation as described in Eqs. (10) and (11) is simplified, such that it can be solved without an iterative procedure. Its initial value equals obs,0 and it monotonically increases with the accumulated plastic strain: Here, material parameters obs and are the hardening factor and exponent. The evolution of the equivalent plastic strain is described as: In Eq. (29), the contribution due to the kink-pair slip resistance follows a similar expression as given in Eq. (12): Here, the parameter values for exponents and are assumed to take the same values as used at the microscopic scale. The formulation regarding the knee-temperature on the macroscopic scale follows a similar expression as given in Eq. (13) as well: Again, the ''knee-strain-rate'' is given bẏe q p,k =̇k p,0 exp(− 0 ∕ B ), and = 10 −7 s −1 . In order to obtain the macroscopic parameter values, the representative volume element (RVE) discussed earlier is uniaxially loaded at various temperatures. Here, the mechanical behaviour is described by the crystal plasticity model given in Section 2.1, with the material parameters listed in Table B.1. In Fig. 4, the blue lines show the homogenized equivalent second Piola-Kirchhoff stress as a function of the equivalent Green-Lagrange strain at various temperatures. The corresponding parameters are obtained in earlier work [7] and summarized in Appendix B.
The macroscopic deformation gradient tensor̄and the macroscopic first Piola-Kirchhoff stress tensor̄are the volume-averages   Fig. 4 shows that with this set of parameters, the macroscopic material behaviour accurately captures the homogenized behaviour of the RVE. For Cu and CuCrZr, isotropic Von Mises plasticity is assumed without hardening and the material parameters in Appendix B.
The increase in yield stress due to neutron irradiation irr is determined for a wide range of irradiation temperatures. This increase in yield stress follows from the increase in slip resistance due to the presence of microscopically small voids. The concentration of these voids follows from the cluster dynamics model as described in Section 2.2. In this analysis, the defect sink strength of grain boundaries is disregarded, as justified in previous work [3].

Failure probability model
In literature, material failure in a monoblock is frequently observed [23][24][25][26], during which a crack propagates in an unstable manner. With Beremin's [22] weakest-link theory, the probability that failure occurs in the monoblock can be assessed. To this end, the monoblock volume is divided into a large number of subvolumes. The probability that failure occurs in the monoblock is governed by the probability that failure occurs in any of its subvolumes: Here, the approximation implies that the probability of failure of each of the subvolumes is small, even when failure in the monoblock is likely to occur.
Following the approach from earlier work [7], the probability of failure of these subvolumes is described by means of a representative volume element (RVE), i.e. f,sub = f,rve . Furthermore, the monoblock geometry is divided into finite elements where the volume of element is denoted as , where = 1, … , elmts . The failure probability of the  monoblock is then expressed as: where rve is the volume of the RVE. In order to describe the probability of failure of the RVE, the weakest link theory is applied once again, and is governed by the probability that failure occurs in any of the grains underlying grains: Here, the approximation is again valid because the failure probability of the RVE is small. The probability of failure of a grain is related to the level of stress subjected to it. According to Griffith's energy balance [2], a crack of length will grow in an unstable manner if it exceeds a critical crack length c = ( Grif f ∕ ⟨ n ⟩) 2 , where Grif f is a constant and ⟨ n ⟩ is the positive normal stress acting on the crack tip opening. Considering a material that is initially flawed by a distribution of microscopically small cracks, the probability of finding a crack of length between and + d within a reference volume ref is assumed to follow a power law ( ) d = ∕ d . Here, and are material parameters. The probability of failure of a grain is governed by the probability of finding a crack that exceeds a critical crack length for a given stress state: Here, 0 is a pre-exponential factor, ref is a measure for the resistance against cleavage and f is a measure for the variation in micro-crack lengths in the material. The reference volume ref is taken equal to the mean grain volume. Cleavage in tungsten is almost exclusively observed to initiate on the {100} family of planes [27]. The assumption that the micro-crack distribution on each of these planes is statistically comparable allows for the following formulation of the failure probability of a grain: Here, the positive grain-averaged Cauchy normal stress can be obtained by means of the following three steps: n , otherwise , n = ⃗ cl ⋅ g ⋅⃗ cl and g = 1 Here, ⃗ cl is the normal vector of cleavage plane , with = 1, 2, 3.
The stress field (⃗) in an RVE is obtained by means of the microscopic material model highlighted in Section 2. Substituting the positive grain-averaged Cauchy normal stress into Eq. (38) yields the failure probability of the grains, which allows for calculation of the failure probability of an RVE with Eq. (36). Finally, the failure probability of the monoblock can be obtained by means of Eq. (35). Note that at the length scale of the monoblock, the notion of ''failure probability'' refers to the probability that a pre-existing crack grows in an unstable manner. This does not necessarily imply that the monoblock can no longer be operational. Furthermore, the crack can get arrested since the geometry and loading conditions that correspond to this boundary value problem are such that the stress intensity at the crack tip does not necessarily increase upon crack propagation. Nevertheless, it will affect the local integrity and heat conduction of the monoblock.

Results
The thermo-mechanical analysis regarding the performance of the monoblock under fusion conditions is divided into two parts. First, the response of the monoblock due to heating is assessed in Section 5.1. Then, two different cooling scenarios are considered in Section 5.2. Here, the impact of irradiation is examined as well.

Heat load
In this analysis, the tungsten monoblock is exposed to a heat load of 10 MW∕m 2 for the duration of 10 seconds, after which the temperature has almost reached a steady state condition. This temperature field is visualized in Fig. 5. The figure reveals a steep temperature gradient in the monoblock.
The temperature distribution in the monoblock induces a spatial variation of the thermal expansion. While areas near and below the cooling tube remain relatively cool, the areas near the plasma facing surface become extremely hot which results in a graded thermal expansion. This mismatch in thermal expansion entails high stresses in the material. Furthermore, the mismatch in the coefficient of thermal expansion between copper and tungsten is another cause of stress near the centre of the monoblock. Finally, the pressure of the water inside the cooling tube also provides a small contribution to the overall stress within the monoblock.
The first principal stress is visualized in Fig. 6 at two different time instances, whereby the highest stress levels are observed already after 1.42 s. The figure reveals that the most severe stress levels are close to the tungsten-copper interface and directly above the cooling tube. Since fracture is mostly observed in the tungsten part of the monoblock, the stress in the cooling tube and interlayer are not explored in the remainder of this work. In the figure, there is no nodal averaging applied. Furthermore, the stresses are plotted by translating the stress values from the integration points to the nodes.
In Fig. 6, the line path A-B is highlighted, and runs along the central line located at the intersection of the two symmetry planes previously indicated in Fig. 3. Point A is located at the plasma facing surface, while point B is located at the W-Cu interface. Along path A-B, crack formation is frequently observed in literature [23][24][25][26] and is therefore examined into closer detail next.
In Fig. 7, the temperature history along path A-B is shown. Here, the -coordinate indicates the depth starting from point A where = 0 mm to point B at which = 8 mm. The figure shows a high temperature gradient, as the temperature varies by about 1000 K over a thickness of only 8 mm. At point A, temperatures above 1300 K are reached. This can trigger recrystallization, which is a phenomenon during which new grains are formed that are nearly defect free. Near point B, the temperature remains relatively low. Tungsten is notorious for its brittleness at low temperatures, where the combination with high stress levels may result in brittle cleavage, which is obviously undesirable as it occurs extremely sudden and can make the monoblock unfit for being operational.
In Fig. 8, the first principal Cauchy stress history along path A-B is shown. The highest stress levels occur close to the cooling tube (i.e. point B), already after 1.42 s, i.e. shortly after the heat source is activated. The temperature profile has not yet reached steady state at this time instant. In this analysis, the focus is on tensile stresses (i.e. the largest eigenvalue of the Cauchy stress tensor), because compressive stresses do not induce material failure.

Microscopic analysis
In order to investigate the response at the microscopic scale, the history of the deformation gradient tensor and the temperature along this path are applied to the representative volume element presented in Section 3.2. Fig. 9 shows the homogenized first principal Cauchy stress evolution that follows from this analysis. As can be seen, it closely reproduces the results obtained from the macroscopic simulation shown in Fig. 8. Note that the results are not expected to be identical due to the simplifications made in the macroscopic model.
The probability of failure of the material is governed by the normal stress acting on the cleavage planes of each grain as described in Eq. (39). In Fig. 10, these stresses are visualized at point B which is located at the W-Cu interface. The results shown in the figure correspond to the stress after 1.42 s of exposure, at which the highest stress levels occur. Since the highest stress levels occur in the width direction (WD) of the monoblock, grains that have cleavage plane normals oriented in this direction endure the most severe loading.

Monoblock failure probability
With the failure probability model discussed in Section 4, the failure probability of the representative volume element (RVE) is calculated, using the Cauchy normal stresses acting on the cleavage planes. To this end, the material parameters 0 = 0.01 and ref = 920 MPa are adopted from earlier work [7], and f = 5.5 is taken from Margetson and Sherwood [28]. The reference volume is again taken equal to the mean grain volume ref = 1.64 × 10 6 μm 3 [7]. The history of the failure probability of the RVE along the path A-B is shown in Fig. 11(a). Clearly, the RVE is most prone to failure at the location and time instant that correspond to the highest stress levels. Fig. 11(a) shows a strong correlation with the homogenized first principal Cauchy stress shown in Figs. 8 and 9. Since the RVE consists of many grains, each with a random crystal orientation, the failure probability model discussed in Section 4 can be closely approximated with the original model from Beremin et al. [22]: With material parameters 0 = 0.01 and f = 5.5 taken equal to those used in the microscopic scale model, a least-squares fit yields ref = 678.0 MPa. With the original model, the failure probability follows directly from the macroscopic first principal Cauchy stress. The model implies that if a macroscopic material volume equal to that of the RVE is subjected to a first principal stress equal tō1 =̄r ef , its failure probability equals 1%. The result of the original model of Beremin et al. shown in Fig. 11(b) reveals a close resemblance with the model results shown in Fig. 11(a).
With the original model from Beremin et al. [22] in Eq. (40), the failure probability history of the entire monoblock is computed by means of Eq. (35), while taking into account that only a quarter of the monoblock geometry is simulated in the finite element calculation. The failure probability during heating of the monoblock is shown by the red line in Fig. 12. As can be seen in the figure, the highest probability of failure in the monoblock is reached before the temperature profile reaches its steady state condition, and can reach a value of around 18%.

Cooling
After a 10 s exposure of 10 MW∕m 2 , the heat load is removed. In Fig. 13, the first principal stress field in the monoblock is shown that corresponds to time = 20 s. As can be seen, the overall stress level is considerably lower compared to the state with a heat load. In Fig. 14(a), the temperature history that corresponds to path A-B is shown. Note that the final temperature profile in Fig. 7 corresponds to the initial temperature in Fig. 14(a). The figure reveals that along this path, the monoblock quickly cools down to the temperature of the cooling fluid.
In Fig. 14(b), the history of the RVE homogenized first principal Cauchy stress is plotted along path A-B. In this figure, the initial stress level corresponds to the final stress level in Fig. 8 and the same colorbar is adopted. The figure reveals that during cooling down, the stress remains relatively low. The stress level does not reduce to zero, because of the plastic deformation that occurs during the heating process. The highest stress level now occurs close to the plasma facing surface (i.e. point A) rather than the W-Cu interface (i.e. point B), which is opposite to what was previously observed in Fig. 8. In other words, the tensile and compressive regions have switched. The failure probability is governed by the stress the material is subjected to. The tradeoff between the brittle failure stress and the temperature-dependent yield stress governs the ductile-to-brittle transition and is captured in this modelling approach. Since the level of stress remains low, the probability of failure of the monoblock remains relatively low as well, which is indicated by the blue line in Fig. 12.

Stress relaxation
In order for nuclear fusion to be economically viable, the fusion process must be operational for a long time period. During this time, recovery processes occur that ultimately result in stress relaxation [29]. In this section, the limit case in which the monoblock is fully relaxed under the static heat load is considered. This is a limit case, as in practice some residual stresses may still remain. Here, the stress relaxation process itself is not simulated. At time relax , the temperature distribution equals the one presented in Fig. 5, but the heat source acting on the stress-free monoblock is turned off and the monoblock rapidly starts to cool down. Under these conditions, the temperature follows virtually the same history as before. However, the magnitude of the stress is considerably higher, as can be seen in Fig. 15 where the stress history after time step relax is shown. Along path A-B, the highest stress levels occur close the plasma facing surface. The stresses acting on the cleavage plane normals at point A are visualized in Fig. 16 and are most significant in the width direction (WD).
The first principal stress in the entire monoblock is revealed in Fig. 17, 10 s after the heat source is shut off. The highest values of the first principal stress are now obtained along the path C-D. Here, the direction of the highest principal stress is parallel to HD, i.e. along the line C-D. This path is also parallel to path A-B, but does not intersect with any of the symmetry planes. In Fig. 18, the blue line shows the corresponding evolution of the failure probability of the monoblock, which is calculated in a similar fashion as presented earlier. The figure reveals that the failure probability now reaches values close to 35%, which is significantly higher than what was obtained during heating. Contradictory to what one might expect, the stress relaxation increases the risk of failure upon shut down of the heat source.

Irradiation
Due to neutron irradiation, the yield stress of tungsten typically increases. This increase is strongly dependent on the temperature at which it is irradiated, as both the source term and the interaction coefficients are temperature dependent. The parameters regarding the interaction coefficients are listed in Table B.2. For the source term, the parameter values c , tot and are obtained by linear interpolation between the temperatures 300 K, 1025 K and 2050 K, following the approach of Mannheim et al. [30]. In the current analysis, an irradiation time of irr = 154.3 h is simulated, which corresponds to a fluence = 5 × 10 19 neutrons∕cm 2 , given a flux of = 3.6 × 10 14 neutrons∕cm 2 [31]. Furthermore, a low dislocation density of d = 10 12 m −2 is assumed since recrystallized tungsten is considered. Defect clusters up to a size of 10 3 are taken into account, while the largest defect size in the source term is c = 50. Regarding the pre-factor of the dispersed barrier hardening model, irr = 0.0314 is taken.
In Fig. 19(a), the increase in the yield stress due to neutron irradiation is visualized. Clearly, neutron irradiation generally increases the yield stress in regions where the temperature is relatively low. In the simulation, the irradiation induced increase in yield stress remains constant once the heat source is shut down, because at this point the neutron source has vanished as well and the temperature throughout the monoblock decreases monotonically.
In Fig. 19(b), the magnitude of the equivalent plastic strain is shown, 10 s after shutting down the heat source. The figure reveals that most of the deformation remains elastic, which implies that the yield stress is not reached. The first principal stress throughout the monoblock is shown in Fig. 19(c). As can be seen, the stress distribution shows a strong resemblance with Fig. 18, although the magnitude of the stress is higher. An even higher fluence hardly yields higher stress levels, since most of the deformation is already elastic. In Fig. 18, the yellow line indicates the failure probability of the irradiated monoblock. Obviously, the failure probability is about 5% higher compared to the unirradiated tungsten that is cooled under the same conditions.

Discussion and conclusions
The process of fusing heavy hydrogen isotopes imposes extreme conditions upon the plasma facing components, among which the highpurity tungsten monoblocks. The large temperature gradient inside  these monoblocks during operation induces large stresses, that may result in a loss of structural integrity.
In order to assess the impact of the loading conditions, a twoscale modelling framework was presented. The tungsten monoblocks have a polycrystalline microstructure that was analysed in detail by means of a representative volume element. At the scale of individual grains, a crystal plasticity model was adopted that describes their mechanical response. Here, the orientation of slip systems is taken into account, leading to an anisotropic stress response. The resistance to crystallographic slip originates from different sources and depends on temperature, the history of plastic deformation, and irradiation history. The latter was described by means of a cluster dynamics model in which the evolution of defect concentrations is governed by a set of rate equations. These equations include the defect production rate, as well as the interaction reactions between defect clusters of various sizes. Voids constitute the most significant obstruction to dislocation propagation and their concentrations pose an additional term in the resistance to crystallographic slip.
The considered tungsten microstructure is fully recrystallized and its texture consists of crystal orientations with a uniformly random distribution. An isotropic Von Mises type of plasticity model was employed that mimics the homogenized response of the polycrystalline microstructure. Although recrystallization at the top layer of the monoblock can be expected already for short durations of high heat input, an interesting topic for further research would be to investigate the performance of rolled tungsten monoblocks. In that case, the texture would be more pronounced and the assumption of isotropy may no longer suffice. In the current work, the evolution of the Von Mises yield stress follows a similar description as the resistance to crystallographic slip. The boundary value problem considered in the analysis of the monoblock consists of quarter of the monoblock geometry, exploiting the two planes of symmetry.
In Beremin's weakest-link theory [22], the failure probability of a given volume is governed by the probability that failure occurs in any of its subvolumes. In this work, this theory was applied on both the macroscopic and microscopic scale. Regarding the latter, the failure probability of the representative volume element is governed by the probability that failure occurs in any of its constituting grains. The probability of failure initiating at a grain boundary is disregarded, as justified in earlier work [7]. However, it remains difficult to experimentally prove that the initial fracture site is always transgranular. The competition between transgranular and intergranular failure should be studied into more detail and could be included in the model in future work. The failure probability of a grain was described by means of Griffith's energy balance [2], and takes into account the orientation of the {100} cleavage planes.
Exposing the water-cooled monoblock to a heat load of 10 MW∕m 2 on the plasma-facing surface results in a significant temperature gradient of more than 1000 K over an armour thickness of 8 mm. A substantial magnitude of stress follows, in particular close to the W-Cu interface. The heat load was applied instantaneously, quickly resulting in high stress levels while the temperature has not yet reached a steady state. This suggests that high stress levels in the monoblock are susceptible to a precipitous temperature profile, whereby the impact of transients and edge localized modes (ELMs) was not even considered in the present analysis. An ELM is a short instability during which several MJ∕m 2 are released.
At several critical points, the macroscopic history of deformation and temperature was applied to a representative volume element. This allowed to analyse the response at the microscopic scale, in which the loading on individual cleavage planes was investigated. The failure probability depends on the presence of initial defects. The material parameters that correspond to this model are obtained from various experiments in literature. In the real monoblock, these parameters could be different due to an incompatible loading history. Controlling the initial defect distribution in the material can be critical in terms of the monoblock performance.
Results showed that the failure probability at the microscopic scale remains consistent with the original model of Beremin et al. [22]. This model is based solely on the macroscopic first principal stress, which can therefore be expected to be insufficient when hydrostatic stress dominates, or when the orientation of the cleavage planes is not randomly distributed.
With the original model from Beremin et al. [22], the failure probability can be assessed at the scale of the monoblock. This refers to the probability that a pre-existing crack or defect initially grows in an unstable manner. For the considered boundary value problem, the stress intensity at the crack tip does not necessarily increase upon crack propagation. Therefore, depending on its location, the crack can get arrested and the monoblock is not necessarily unfit for operation. Shortly after exposure to the heat load, the monoblock failure probability rises to about 18%, after which it monotonically decreases to about 8%. These are particularly high values, given the fact that the fusion reactor contains hundreds of thousands of monoblocks. The highest magnitude of stress is situated close to the W-Cu interface, which is therefore also the site at which failure is most likely to initiate.
Once the heat load is removed, the monoblock quickly cools down to the reference temperature. Furthermore, the stress magnitude drops dramatically, but does not reach zero because of the plastic deformation that occurred during prior heating. In contrast to what was observed during heating, the highest stress levels are now located close to the plasma facing surface. The overall failure probability of the monoblock during cooling remains relatively low. However, this perception drastically changes if the stress relaxation inside the monoblock is assumed to be complete under the applied static heat load. Now, the failure probability reaches values close to 35%, but the highest stress levels occur at the corner edges of the monoblock. Fig. 18 has shown that within the current modelling framework, the impact of neutron irradiation on the failure probability is limited to 5%. This follows from the fact that most deformation is already elastic. The simulation results show that plastic deformation occurs solely close to the plasma facing surface, where the impact of neutron irradiation is minimal due to the high irradiation temperature that corresponds to this location. In this analysis, the decrease in thermal conductivity due to neutron irradiation is not taken into account. This aspect is expected to increase the temperature in the monoblock, which also yields an elevated amount of thermal expansion and hence a higher stress level. Furthermore, the impact of neutron irradiation on the yield stress of copper and the CuCrZr cooling tube is also not taken into account. The yield stress is known to increase upon irradiation, and it is therefore expected that increased stress levels will emerge close to the W-Cu interface. In future research, it would be interesting to incorporate this aspect by including an increased yield stress that corresponds to the same irradiation conditions considered here.
Experimental data in literature suggests that for tungsten, the resistance to cleavage increases upon neutron irradiation [32,33]. Phenomena such as dislocation channelling suggest that the magnitude of plastic deformation becomes increasingly heterogeneous after irradiation [34,35]. As pointed out in previous work [3], the mechanism causing an increase in cleavage stress due to neutron irradiation is not yet resolved. Here, the brittle-to-ductile transition temperature of irradiated specimens is described by assuming an enhanced cleavage resistance of ref = 1.8 GPa, compared to ref = 920 MPa for the unirradiated specimens. Such a drastic increase in the cleavage resistance can be expected to drastically reduce the failure probability. In this regard, future research should be directed towards the understanding of the underlying mechanism that governs this behaviour.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix B. Material parameters
The material parameters required for the crystal plasticity model in Section 2.1 and the cluster dynamics model 2.2 are listed in Tables B.1 and B.2, respectively.
Temperature dependent Young's modulus, specific heat, thermal conductivity and yield stress are adopted from the ITER materials handbook [45], and visualized in Fig. B.20. Furthermore, the mass densities of copper and CuCrZr are also adopted from this reference. These are taken to be 8.940 g∕cm 3 and 8.900 g∕cm 3 , respectively. The thermal expansion coefficient is taken to be W = 3.934 × 10 −6 K −1 and Cu = CuCrZr = 16.84 × 10 −6 K −1 [45]. Finally, the Poisson ratio of tungsten is taken to be = 0.29, and = 0.33 for Cu and CuCrZr, respectively [45]. Fig. B.20. Temperature dependent material parameters adopted from the ITER materials handbook [45]. The yield stress of tungsten is not only dependent on temperature, but also on the equivalent plastic strain and plastic strain rate, as highlighted in Section 3.2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Table B.1 Material parameters of tungsten for the crystal plasticity model, adopted from Asaro [8], Brunner [36], Kocks et al. [37], Lim et al. [38], Lowrie and Gonas [39], Yalçinkaya et al. [13] and Oude Vrielink et al. [3,7,9].

Parameter
Symbol Value