Yixiang Gan

Francisco Hernandez

Dorian Hanaor

Ratna Annabattula

Marc Kamlah

Pavel Pereslavtsev

A THERMAL DISCRETE ELEMENT ANALYSIS OF EU SOLID BREEDER BLANKET SUBJECTED TO NEUTRON IRRADIATION

Keywords: tritium breeding materials, test blanket module, discrete element method, thermo-mechanical analysis, neutron heating, granular material

Due to neutron irradiation, solid breeder blankets are subjected to complex thermo-mechanical conditions. Within one breeder unit, the ceramic breeder bed is composed of sphericalshaped lithium orthosilicate pebbles, and as a type of granular material, it exhibits strong coupling between temperature and stress fields. In this paper, we study these thermo-mechanical problems by developing a thermal discrete element method (Thermal-DEM). This proposed simulation tool models each individual ceramic pebble as one element and considers grainscale thermo-mechanical interactions between elements. A small section of solid breeder pebble bed in HCPB is modelled using thousands of individual pebbles and subjected to volumetric heating profiles calculated from neutronics under ITER-relevant conditions. We consider heat transfer at the grain-scale between pebbles through both solid-to-solid contacts and the interstitial gas phase, and we calculate stresses arising from thermal expansion of pebbles. The overall effective conductivity of the bed depends on the resulting compressive stress state during the neutronic heating. The thermal-DEM method proposed in this study provides the access to the grainscale information, which is beneficial for HCPB design and breeder material optimization, and a better understanding of overall thermo-mechanical responses of the breeder units under fusion-relevant conditions. 1

I. INTRODUCTION

During operations of nuclear fusion breeder reactors, Helium-Cooled Pebble Beds (HCPB), as a candidate for fusion blankets, [START_REF] Boccaccini | The EU TBM systems: Design and development programme[END_REF] are subjected to neutron fluxes to generate tritium to complete the fuel cycle. Thermo-mechanical conditions inside the breeder units are crucial for tritium generation and release rate, [START_REF] Ying | Status of ceramic breeder pebble bed thermo-mechanics R&D and impact on breeder material mechanical strength[END_REF] and the mechanical stability [START_REF] Gan | Crush probability analysis of ceramic breeder pebble beds under mechanical stresses[END_REF] of the ceramic breeder materials, e.g., lithium orthosilicate, Li 4 SiO 4 . [START_REF] Knitter | Crystallisation and microstructure of lithium orthosilicate pebbles[END_REF] In order to gain an accurate understanding of the conditions during operations, thermo-mechanical models and experiments have been carried out for those materials, [START_REF] Reimann | Thermomechanical behaviour of ceramic breeder and beryllium pebble beds[END_REF][START_REF] Reimann | Thermomechanics of solid breeder and Be pebble bed materials[END_REF][START_REF] Reimann | Thermal conductivity of compressed ceramic breeder pebble beds[END_REF][START_REF] Dell'orco | Progress in the benchmark exercise for analyzing the lithiate breeder pebble bed thermo-mechanical behaviour[END_REF][START_REF] Zhao | Influence of plate material on the contact strength of Li4SiO4 pebbles in crush tests and evaluation of the contact strength in pebbleâ Ȃ Şpebble contact[END_REF] including phenomenological and grain-scale approaches.

Due to its granular nature, the solid breeder material exhibits strong thermo-mechanical coupling and heterogeneous behaviour. Experimental data [START_REF] Reimann | Thermomechanics of solid breeder and Be pebble bed materials[END_REF][START_REF] Reimann | Thermal conductivity of compressed ceramic breeder pebble beds[END_REF][START_REF] Abou-Sena | Effective thermal conductivity of lithium ceramic pebble beds for fusion blankets: A review[END_REF] show a clear straindependent thermal conductivity, which raises questions concerning both the understanding of the underlying physics and [START_REF] Boccaccini | The EU TBM systems: Design and development programme[END_REF] Fusion Science and Technology, DOI: 10.13182/FST13-727. the implementation of quantitative engineering analysis by considering this coupled thermo-mechanical behaviour. In phenomenological models, this coupling behaviour has been implemented based on empirical curves for effective thermal conductivity of ceramic beds from elaborated experimental measurements. The experimentally obtained effective thermal conductivity can be expressed as a function of both bed temperature (T ) and strain (ε), i.e., k eff (T, ε). [START_REF] Reimann | Thermomechanics of solid breeder and Be pebble bed materials[END_REF][START_REF] Reimann | Thermal conductivity of compressed ceramic breeder pebble beds[END_REF] However, the empirical curves and coefficients may vary if different breeder materials, for instance lithium metatitanate, will be used or simply adopting different types of pebble size distributions of pebbles. Without accessing grain-scale information under operation-relevant conditions, the thermo-mechanical behaviour of the material can not be accurately assessed within breeder units.

Combining discrete and continuum approaches preserves both micro-scale material details and modelling capability at the system scale. Achieving this requires a multi-scale modelling scheme for bridging information at two different length scales. Statistical data can be achieved under given thermo-mechanical conditions, e.g., the force acting on individual pebbles, [START_REF] An | Numerical characterization of thermo-mechanical performance of breeder pebble beds[END_REF][START_REF] Gan | Discrete element modelling of pebble beds: With application to uniaxial compression tests of ceramic breeder pebble beds[END_REF][START_REF] Annabattula | Mechanics of binary and polydisperse spherical pebble assembly[END_REF] and the crushing of pebbles. [START_REF] Gan | Crush probability analysis of ceramic breeder pebble beds under mechanical stresses[END_REF][START_REF] Zhao | Influence of plate material on the contact strength of Li4SiO4 pebbles in crush tests and evaluation of the contact strength in pebbleâ Ȃ Şpebble contact[END_REF][START_REF] Annabattula | Mechanics of a crushable pebble assembly using discrete element method[END_REF] To understand the thermo-mechanical coupling in the material with a granular form, we propose a numerical method, i.e., thermal discrete element method (Thermal-DEM [START_REF] Gan | Discrete element modelling of pebble beds: With application to uniaxial compression tests of ceramic breeder pebble beds[END_REF][START_REF] Gan | Phase transitions and cyclic pseudotachylyte formation in simulated faults[END_REF] ), to estimate the temperature and stress profiles within the ceramic breeder material under ITER-relevant conditions. Thermal-DEM models ceramic breeder pebbles as individual elements and describes the interactions between neighbouring pebbles, including not only contact forces but also heat fluxes, via solidto-solid contacts and the interstitial gas phase between pebbles. A typical EU design for breeder units [START_REF] Hernández | Set-up of a pre-test mock-up experiment in preparation for the HCPB Breeder Unit mock-up experimental campaign[END_REF] and the relevant neutronic heating profile [START_REF] Pereslavtsev | Neutronic analysis of the HCPB TBM in ITER utilizing an advanced integral approach[END_REF] have been adopted for the numerical analysis to demonstrate the feasibility of the proposed method for the analysis of solid breeder blanket systems. By identifying key parameters that determine the resulting temperature profiles and stresses, this method can be incorporated for the design and optimisation of breeder materials in fusion blankets.

II. METHOD

II.A. Thermal Discrete Element Method

To study the heat transfer in pebble bed systems, we use the thermal discrete element method (Thermal-DEM), which is based on the conventional DEM [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] with one additional degree of freedom for the temperatures of the pebbles, to model individual pebbles and their thermo-mechanical interactions.

The mechanical interactions between pebbles, which are modelled as elastic spheres, include the normal and frictional forces, based on the analytical solution for Hertzian contact. The positions of pebbles are updated according to Newton's law of motion with an explicit integration scheme, see more details in. [START_REF] Gan | Discrete element modelling of pebble beds: With application to uniaxial compression tests of ceramic breeder pebble beds[END_REF] For updating the temperatures of individual pebbles, the heat conduction between pebbles is described through Fourier's law. The thermal conductance between the pebbles consists of two parts: solid contact conductance, and the gap conductance through the interstitial gas phase. The inter-granular conduction is described by the analytical model proposed by, [START_REF] Batchelor | Thermal or electrical conduction through a granular material[END_REF] including the combined modes of heat transfer through both contact regions and gap regions. This analytical model considers the relative thermal conductivity between the solid and gas phases, as k s /k g , where k s denotes the conductivity of the bulk material and k g is the gas conductivity. Moreover, the heat generation inside each pebble can be calculated by the power spectrum based on the neutron irradiation profile. This profile is introduced in Thermal-DEM as a source of volumetric heating.

Combining the heat transfer and heat generation by neutron irradiation, the rate of temperature change ( Ṫi ) for the i-th pebble can be updated as Ṫi = 1

m i c p ( j φ i j + Ψ * V cell i ) . (1)
Here, m i and c p are the mass of the pebble and the heat capacity of the bulk material, respectively. The heat flux φ i j is calculated by

φ i j = -2πk g (T i -T j )R eff H(α, β) . ( 2
)
The non-dimensional flux H is based on the relative thermal conductivity, α = k s /k g , and relative contact area, β = αa/(2R eff ), where a is the radius of contact area and

R eff = R i R j (R i + R j ) -1
is the effective radius (more details can be found in [START_REF] Batchelor | Thermal or electrical conduction through a granular material[END_REF] ).

We adopt the volumetric heating profile from the neutronics calculation, [START_REF] Pereslavtsev | Neutronic analysis of the HCPB TBM in ITER utilizing an advanced integral approach[END_REF] as Ψ * shown in Fig. 1(b). The actual heating power in each pebble can be scaled based on its local volume and the power density, as Ψ * V cell i , where the cell volume, V cell i , is obtained using a Voronoi tessellation for individual pebbles, [START_REF] Rycroft | Analysis of granular flow in a pebble-bed nuclear reactor[END_REF] shown schematically as the dashed cell in Fig. 1(c). For each pebble, the scaled heating power density per solid volume is given by

Ψ i = Ψ * V cell i /V pebble i
, where

V pebble i
is the solid volume of the i-th pebble.

Thermo-mechanically coupled behaviours in the grainscale model are automatically taken into account through the thermal expansion of pebbles. Under conditions of neutron heating, a compressive stress develops in the bed as the result of the pebbles' thermal expansion. In turn, the growing compressive stress enhances heat transfer through the enlargement of contact areas.

Based on the EU breeder unit design [START_REF] Hernández | Set-up of a pre-test mock-up experiment in preparation for the HCPB Breeder Unit mock-up experimental campaign[END_REF] in Fig. 1(a), the neutron heating profile for the ceramic breeder materials has been calculated in [START_REF] Pereslavtsev | Neutronic analysis of the HCPB TBM in ITER utilizing an advanced integral approach[END_REF] for Li 4 SiO 4 . To demonstrate the feasibility of using this numerical method, three power densities have been selected, namely, 5, 7 and 8 MW/m 3 , corresponding to different radial positions inside the breeder layer. The exact positions are not given in this study, since only the heating power density will be required from the blanket design.

II.B. Material Parameters and Boundary Conditions

The material properties of Li 4 SiO 4 pebbles, compiled from different sources, are listed in Table I. Those values have been calculated only at a fixed temperature of 500 • C due to their weak temperature dependencies within the targeted temperature range. It would be straightforward to implement bulk material properties varying with temperature in the future work.

The interstitial gas properties are listed in Table II. Note that, relative to the bulk material, the thermal conductivity of the gas phase is sensitive to temperature. In this study, temperaturedependent gas conductivity has been implemented according to the empirical equation provided by. [START_REF] Petersen | The properties of helium: density, specific heats, viscosity, and thermal conductivity at pressures from 1 to 100 bar and from room temperature to about 1800 K[END_REF] The empirical equation used in this study includes also pressure dependency, but the influence of gas pressure variation from 100 to 1000 kPa is negligible. A fixed purged gas pressure of 400 kPa is used in our simulations. It is noted that from experimental observations a higher interstitial gas pressure can lead to better bed conduction. [START_REF] Abou-Sena | Effective thermal conductivity of lithium ceramic pebble beds for fusion blankets: A review[END_REF] This is mainly due to the heat conduction via the gas phase in the regime of large Knudsen numbers (K n = λ/L) depends on both the dimensions of the gaps (L) and the mean free path of the gas molecules (λ). In a granular media, the gap size (L) can be linked to the size of pebbles. This effect is not considered in this study, but it clearly will influence the effective conductivity of the beds, especially for beds with small pebbles.

The height of the simulation cell is 22 mm (along the xaxis, adopted from the EU blanket design shown in Fig. 1(a)), and its width is 5 mm along the y-axis. The depth (the z-axis) is varied to accommodate 5,000 pebbles at a given initial packing factor, η, but the depth dimension is more than 5 times larger than the mean diameter of pebbles. For y and z directions, periodic boundary conditions are applied, while at the x direction elastic wall conditions have been implemented with properties of EUROFER listed in Table II. The wall has a constant temperature of 500 • C during the neutron heating stage. Prior to the neutron heating, a uniform heating stage is performed to raise the system temperature from 20 • C to 500 • C for all the pebbles and walls simultaneously. Then, the neutron heating stage starts with a constant heating profile Ψ * . Two kinds of size distributions have been considered inside the simulation cells: (1) mono-sized samples with diameters of 0.4±0.02 mm; (2) polydisperse samples with diameters ranging from 0.25 to 0.65 mm, with a normal distribution with respect to pebble volume. Samples with various initial packing factors have been generated at 20 • C with a negligible level of compressive stress, i.e., p = -σ ii /3 < 0.1 MPa.

To transfer the physical properties of materials listed in Table I and Table II into DEM simulations, some simple scaling laws have been applied to speed up the calculation. Three time scales are considered here: (1) collision time, t c = ρd 2 /E, describes the typical time for collision events between two elastic pebbles, with elastic modulus of E, diameter of d and density of ρ; (2) thermal diffusion time, t th = ρc p d 2 /k s , describes the characteristic time to transfer the heat between pebbles, with conductivity of k s and heat capacity of c p ; and (3) heating time, t h = ρc p ∆T/Ψ, provides an estimation for the time to heat a pebble with a temperature increase of ∆T via a given heating power density Ψ. In the current type of problems, the collision time (t c ) determines the time step required to achieve meaningful solutions under the explicit scheme of DEM.

To achieve a larger time step in simulations, a combination of scaled density, size and elastic modulus can be selected by using the scaling laws. Here, the scaling principle is to ensure the same ratio between the thermal time t th and heating time t h for both physical and simulation domains. Therefore, in terms of thermal analysis, the time used in the simulation corresponds to a unit of time in a physical system. Meanwhile, the condition of t c t th and t c t h ensures the preservation of a mechanical quasi-static state during the heat transfer and heat generation. More details can be found in. [START_REF] Gan | Phase transitions and cyclic pseudotachylyte formation in simulated faults[END_REF] The conversion between physical units and simulation units can be found in Table III.

III. RESULTS

In our DEM simulations, samples with two kinds of size distributions of pebbles are considered: (1) mono-sized beds with an initial a packing factor ranging from 62.5% to 64.5%;

(2) polydisperse beds with packing factor from 63.5% to 66.0%. After the initial uniform heating from 20 • C to 500 • C, three different heating power densities, Ψ * =5, 7, and 8 MW/m 3 , are applied uniformly across the simulation cell until the cell reaches the steady state. In such cases the variation of the power density along radial direction has been ignored. This can be justified since the ratio between the width and height of the cell is sufficient small, and thus the temperature variation along width is less significant compared with the variation along the cell height. However, under this simplification, the compressive stresses may be overestimated for regions having higher power density due to the restriction of pebble movement along the radial direction to the neighbouring regions having less thermal stresses.

Fig. 2 shows typical evolutions of the average bed temperature and hydrostatic pressure (p = -σ ii /3) over operation time from the start-up of neutron heating. The steady states are reached after a few hundred seconds after the start-up. Here, only a mono-sized sample with packing factor η=63.5% is shown. Other samples have similar behaviour but with different magnitudes of bed temperature and pressure at the steady state. Prior to the neutron heating stage, the uniform heating from room temperature to 500 • C introduces a hydrostatic pressure of 4.5 MPa inside the bed for this mono-sized sample (η=63.5%), and the pressure develops to different levels at the steady state, depending on the corresponding magnitude of the neutron heating power density. In this analysis, plastic and creep deformations of pebbles have not been considered. The resulting thermal stresses can be reduced if inelastic behaviour is taken into account at the microscopic scale.

Next, the bed temperature and hydrostatic pressure at the steady state (here t=300 sec is selected) for different samples are categorised based on their size distribution of pebbles and plotted against the initial packing factor in Fig. 3. For the same type of size distribution, samples with higher packing factors result in lower bed temperatures but higher stress levels, subjected to a given heating power density. Relative to a mono-sized sample, a polydisperse sample having a similar packing factor can have a lower stress level and a higher bed temperature, under the same neutron heating profile. The stress level can be connected to the forces acting on pebbles. Previous studies provide some scaling laws relating to these two parameters. [START_REF] Gan | Discrete element modelling of pebble beds: With application to uniaxial compression tests of ceramic breeder pebble beds[END_REF][START_REF] Annabattula | Mechanics of binary and polydisperse spherical pebble assembly[END_REF] The correlations between resulting temperature, stress and packing factor provides information for the selection of materials and the design of components for various loading conditions. For instance, the resulting stress state is highly relevant for the determination of the quantitative requirements for crush loads of pebble materials. The average and maximum temperature can be estimated as a key parameter for blanket design and optimisation for tritium breeding and release. In these cases, the grain-scale information provided by DEM proposed here is of great value. Fig. 4 shows the temperature distribution at the steady state for different size distributions and initial packing factors. In Fig. 4, the volume averaged temperature profiles are shown with respect to vertical position, across the ceramic pebble layer. With the same initial packing factor, for example η=0.635 and 0.640, the polydisperse sample presents lower effective conductivity and lower compressive stress level, and thus a higher bed temperature, as compared to the corresponding mono-sized sample, as shown also in Fig. 3(c). Some "hot spots" are observed in the temperature distributions, and this may result from two factors: (1) a pebble has almost no solid-to-solid contact with its neighbouring pebbles, which leads to a high thermal resistance for heat to be transferred to its surrounding pebbles;

(2) the heating power is calculated based on Ψ * V cell i , while for some small pebbles in polydisperse samples, in particular in Fig. 4(b), the volume occupied by the solid, V pebble i , is comparably small in relation to the total Voronoi volume V cell i giving rise to high heating power density in those pebbles. The second factor for the occurrence of "hot spots" can be addressed if the neutronics calculation can consider the heterogeneity within the granular materials, and provide a better way to take into account the transition between the average heating profile and the actual heating power for individual pebbles. The maximum temperature found inside the bed is around 650 • C, which is lower compared with the finite element (FE) calculations, [START_REF] Hernández | Set-up of a pre-test mock-up experiment in preparation for the HCPB Breeder Unit mock-up experimental campaign[END_REF] around 720 • C. The main reasons are: (1) the FE analyses have not considered the effective thermal conductivity with strain dependency which overestimated the maximum temperature;

(2) DEM simulations have a pre-heating stage, which results in a hydrostatic pressure of a few MPa and thus a higher bed conductivity even before the neutron heating; and (3) ignoring the region of a large Knudsen Number may result in overestimation of the effective conductivity of the bed. In general, mono-sized and polydisperse samples have different temperature profiles along the vertical direction, depending on the initial packing factor. For polydisperse samples, temperatures near the wall regions are more sensitive to the initial packing factor. Under the same hydrostatic stress level and heating profile, mono-sized samples have larger temperature gradients, indicating that the bed has a lower effective thermal conductivity, relative to the polydisperse pebble beds.

IV. CONCLUSION

In this study, a thermal discrete element method (Thermal-DEM) has been developed for ceramic breeder materials under neutron irradiation, which includes volumetric heating, thermomechanical coupling, and inter-granular thermal conduction. The analysis of ceramic breeder materials subjected to neutron heating under ITER-relevant conditions has demonstrated that the resulting stress state and temperature distribution depend on the type of pebble size distribution of the pebbles and their initial packing factor. At the same heating profile and packing factor, as compared to a mono-sized bed, a polydisperse bed can reduce the stress level, which lowers the risk of crushing pebbles in the bed, but it results in a higher bed temperature. Moreover, this study provides a numerical method to estimate thermo-mechanical responses of beds with various packing factors and pebble size distributions, which can be used as an optimisation tool for the design of solid breeder units.

Fig. 1 .

 1 Fig. 1. Schematics of the approach: (a) EU design for ITER-TBM unit adopted from; 16 (b) calculated density of neutron irradiation during operation from; 17 (c) Thermal-DEM model for two pebbles, with mechanical interaction F i j , volumetric heating Ψ i , and heat fluxes φ cont i j and φ gap i j for solid-tosolid and gap conduction, respectively.

Fig. 2 . 3 Fig. 3 .

 233 Fig.2. Evolution of averaged bed temperature (left) and hydrostatic pressure (right) for different heating power densities as a function of neutron heating duration. Sample: mono-sized, η=63.5%.

(Fig. 4 .

 4 Fig. 4. Predicted temperature profiles at the steady state (t=300 sec, ⇤ =8 MW/m 3 ): (a) mono-sized samples, (b) polydisperse samples, with di↵erent initial packing factors. Pebbles are coloured by their corresponding temperature. In the last column, the volume averaged temperature profiles for each of those samples are plotted against the vertical distance (the walls are located at 0 and 22 mm).

Fig. 4 .

 4 Fig. 4. Predicted temperature profiles at the steady (t=300 sec, Ψ * =8 MW/m 3 ): (a) mono-sized samples, (b) polydisperse samples, with different initial packing factors. Pebbles are coloured by their corresponding temperature. In the last column, the volume averaged temperature profiles for each of those samples are plotted against the vertical distance (the walls are located at 0 and 22 mm).

TABLE I .

 I Properties for Li 4 SiO 4 at 500 • C and porosity 5%.

	Property Young's modulus (GPa) 22 Poisson's ratio thermal conductivity (W/mK) 10 heat capacity (J/kg K) 23 thermal expansion coefficient 23 density (kg/m 3 ) pebble size (mm)* pebble-pebble friction pebble-wall friction	Symbol Value E 92 ν 0.24 k s 2.743 c p 2025 α 2.28×10 -5 ρ 2270 d 0.25-0.65 µ 0.2 µ w 0.2
	* pebble size distributions used in simulations include (1) mono-sized samples with 0.4 ± 0.02 mm; (2) polydisperse samples of pebbles between 0.25 and 0.65 mm.

TABLE II .

 II Interstitial gas properties and material properties of the container wall, at 500 • C.

	Property	Symbol Value
	purge gas (helium): thermal conductivity (W/mK) 21 gas pressure (kPa) wall (EUROFER 24 ): Young's modulus (GPa) Poisson's ratio thermal conductivity (W/mK)	k g p g E w ν w k w	0.3 400 175 0.3 29.0

TABLE III .

 III Unit system conversion between physical quantities and simulation parameters.

	Property Young's modulus force pebble size conductivity heat capacity density heating profile temperature	Symbol E F d k s , k g c p ρ Ψ ∆T	Unit 10 9 Pa 10 N 10 -3 m 1 W/mK 1000 J/kg K 1 Simulation 100 1 1 1 ** 1000 kg/m 3 1 10 6 W/m 3 1 ** 1 • C 1
	** calculated via the scaling law between physical system and simulations: t th /t sim th = t h /t sim h .

ACKNOWLEDGMENTS

Financial support for this research from the Australian Research Council through Grant No. DE130101639 is gratefully appreciated. The Thermal-DEM code and scripts for the present study can be downloaded from Github through https://github.com/ganyx/ThermalDEM-TBM.